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The nature of the zero temperature ordering transition in the three dimensional Gaussian random 
field Ising magnet is studied numerically, aided by scaling analyses. Various numerical calculations 
are used to consistently infer the location of the transition to a high precision. A variety of bound- 
ary conditions are imposed on large samples to study the order of the transition and the number of 
states in the large volume limit. In the ferromagnetic phase, where the domain walls have fractal 
dimension da = 2, the scaling of the roughness of the domain walls, w ^ L'^ , is consistent with the 
theoretical prediction (^ = 2/3. As the randomness is increased through the transition, the probabil- 
ity distribution of the interfacial tension of domain walls scales in a manner that is clearly consistent 
with a single second order transition. At the critical point, the fractal dimensions of domain walls 
and the fractal dimension of the outer surface of spin clusters are investigated: there are at least two 
distinct physically important fractal dimensions that describe domain walls. These dimensions are 
argued to be related by scaling to combinations of the energy scaling exponent, 9, which determines 
the violation of hyperscaling, the correlation length exponent i/, and the magnetization exponent /3. 
The value /3 = 0.017 ± 0.005 computed from finite size scaling of the magnetization is very nearly 
zero: this estimate is supported by the study of the spin cluster size distribution at criticality. The 
variation of configurations in the interior of a sample with boundary conditions is consistent with 
the hypothesis that there is a single transition separating the disordered phase with one ground 
state from the ordered phase with two ground states. The array of results, including values for 
several exponents, are shown to be consistent with a scaling picture and a geometric description 
of the influence of boundary conditions on the spins. The details of the algorithm used and its 
implementation are also described. 
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I. INTRODUCTION 



In spite of many years of study, the behavior of phases 
and phase transitions that arc dominated by quenched 
randomness is still controversial. One such lively contro- 
versy has concerned the existence or lack thereof of an 
ordered phase in the random field Ising model (RFIM) in 
three dimensions. Although this waS|-eventually resolved 
in the affirmative by rigorous work,tl the nature of the 
phase transition and the possibility of a phase intermedi- 
ate between the paramagnet and the ferromagnet is still 
controversial. 

Numerical simulations of the random field Ising model 
— and experiments — are impeded by the dramatic 
slowing down that occurs as the phase transition is ap- 
proached due to the existence of free energy barriers 
which are broadly distributed but typically grow as a 
power of the correlation length. Such barriers are gen- 
eral characteristics of phases controlled, in a renormaliza- 
tion group (RG) sense, by stable zero temperature fixed 
points. For the random field Ising model, not only is the 
low temperature phase controlled by a zero-temperature 
fixed point (as is the case for conventional pure systems), 
but the phase, transition itself is also controlled by such a 
fixed pointfl'u Indeed, the ground state properties of the 
random field Ising model undergo a phase transition as 
the strength of the randomness is increased and it is this 



zero-temperature transition that governs the behavior of 
the transition at positive temperatures. Fortunately, this 
means that much can be learned by studying the ground 
state properties. 

It has been known for some time that combinatorial al- 
gorithms can be used effectively to find the ground states 
of various classes of random systems and thft RFIM was 
one of the first to be studied in this wayel With cur- 
rent computers, the algorithm is very fast and large sys- 
tem sizes can be studied in enough detail to obtain good 
statistics enabling the tools of finite size scaling to be 
used to analyze the zero-temperature phase transition. 

Various significant open questions exist about the 
phase transition in the RFIM. Although a self-consistent 
scaling picture of a zero-temperature critical fixed point 
was proposed early on, it has not been adequately tested 
and other scenarios havci-been suggested, including a 
first order phase transitionElia and an. intermediate phase 
with "replica-symmetry breaking"l3Q presumably mean- 
ing many coexisting equilibrium states. 

In this paper we study the RFIM with Gaussian dis- 
tributed random fields, focusing on the nature of the 
phase transition and the sensitivity of the ground states 
to varying boundary conditions (BCs) as a probe of the 
number and nature of the infinite system states. As will 
be explained in some detail, our results strongly support 
the scaling picture of the transition. In this picture, there 
is a single second order critical point characterized by 



three scaling exponents: v for deviations from the criti- 
cal point, Q for the energy at the critical point, and ^ jv 
for the magnetization at the critical point. We clarify 
some of the substantial confusion about the order of the 
transition by showing that both d and v as well as the 
distributions of the "stiffness" and spin clusters are very 
different from what they would be at a first order tran- 
sition. Nevertheless, as observed previously, the magne- 
tization exponent /3 is extremely small so that even with 
the very large sizes we study, the magnetization appears 
almost discontinuous. 



II. MODEL AND NUMERICAL METHOD 

The random field Ising modelB has Hamiltonian, de- 
fined over spin configurations {s^ = ±1}, 
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Y^K 



(1) 



with the random fields hi chosen independently from a 
distribution which we take to be Gaussian with mean zero 
and variance Y? . The ferromagnetic exchange coupling, 
J, is fixed at unity in the simulations and the sites i lie on 
a cubic lattice with interactions between nearest neigh- 
bor pairs {(ij)}. The basic nature of the phase diagram 
of the 3-D RFIM is well-known: As the temperature is 
lowered for small /i, there will be a critical temperature, 
Tc{h), below which the RFIM becomes ferromagnetically 
ordered with a non-zero spontaneous magnetization. As 
the strength of the random field increases, the critical 
temperature decreases until at a critical field, h^, it goes 
to zero. 

Both the paramagnetic and the ferromagnetic phases 
have been ptpven to exist at both zero and positive 
temperatureil and the transition between them can thus 
be studied by varying ft, at T = 0. The simplest sce- 
nario at zero-temperature is a single critical field strength 
he above which the spins are disordered with a unique 
infinite-system ground state and exponential decay of 
correlations, and below which there are two infinite- 
system ground states, one with predominantly up spins 
and the other with predominantly down spins. 

The nature of the phases and the phase transition(s) 
between them can be probed by studying the effects of 
various boundary conditions on larger and larger systems 
— most simply cubes of size V ^ Lx Lx L. In the disor- 
dered phase the orientation of a spin far from the bound- 
aries is typically determined by the collection of random 
fields within a correlation length ^(/i) of the spin and 
is insensitive to boundary conditions imposed far away. 
In contrast, in the ferromagnetic phase some spins will 
still be controlled by the random fields in their vicinity, 
but a finite fraction of the spins will be controlled by the 
boundary conditions — no matter how far away they are 
imposed. The simplest scenario is a single transition be- 
tween these two phases. The primary goal of this paper is 



to examine in detail the nature of this zero-temperature 
phase transition. 

Many previous studies of the ground states of the 
RFIM (as well as finite temperature studies) have focused 
on the magnetization per spin, m = V~^ ^^ s^, and the 
results have been somewhat ambiguous. Some have in- 
terpreted the muTiPTjical results as indicating a second or- 
der transitionJJIlB'EJ whole others have concluded tbat the 
transition is first order .Q Some Monte Carlo resultsti sug- 
gest that the finite-temperature transition is second or- 
der, buj_svith the magnetization exponent f3 nearly zero. 
Othersjlj using a varying external field, have found a co- 
existence of states suggestive of a first order transition. It 
is clear from these studies that if the transition is second 
order, the order parameter exponent /3, \m\ ^ {he — h)P 
must be very small, making definitive conclusions based 
on magnetization alone difficult. We have thus focused 
much of our attention — particularly for locating the 
transition and finding the exponents — on other proper- 
ties which naturally distinguish the phases. 



A. Algorithm 

The (almost surely) unique ground state of a fi- 
nite sample can be datermined in time polynomial in 
the number of spins.Lj The method is based on a re- 
duction of the problem of determining RFIM ground 
states to a maximum-flow problem on an augmented 
graph. O pa caB i t hen use combinatorial optimization 
algorithmalaOOH to solve the maximum-flow problem. 
We describe the special features of the algorithm imple- 
mentation, its verification, sample timings, and the use 
of integer valued hi in the Appendix. 



B. Statistics and analysis 

We have studied system sizes up to 256"^, which contain 
over 1.6 x 10^ spins. Independent samples were simulated 
for each value of h. Separate realizations were also gen- 
erated for boundary induced domain walls, spin cluster 
properties, magnetization, and the thermodynamic limit 
studies. The same samples and domain walls were used 
in the stiffness and domain wall property studies. For 
smaller systems (8^ through 32'^), 10^ samples were op- 
timized, typically. (For the domain roughness measure- 
ments in the ordered phase, 10^ or 10^ samples provided 
sufficient data, as fluctuations in the interface width are 
not large.) Of order 10"^ to 10* samples were studied for 
each quantity for the 64'^ and 128"^ samples. For L = 256, 
4 X 10^ to 10^ samples were studied at each h, as part of 
the magnetization and cluster studies. 

Error bars for exponent values throughout this paper 
include both estimated systematic errors due to appar- 
ent finite size effects and errors due to statistical uncer- 
tainties; the error bars represent an estimated range of 
values in which the value lies, with high confidence. In 



contrast, error bars in the figures reflect la statistical un- 
certainties, which we find to be generally consistent with 
confidence intervals found by resampling. 

Generally (except for the stiffness, the roughness in 
the ferromagnetic phase, fitting a power law to the bond 
energy density, and the Pd/o,± plots), we have used esti- 
mates of effective exponents as a function of system size 
to estimate exponents, rather than scaling plots. This is 
done to more clearly see trends in the data that reflect 
finite size corrections. Finite size corrections tend to be 
monotonic and introduce a drift with L in the effective 
exponents. Given the good statistics of the data sets that 
can be generated with optimization algorithms, collaps- 
ing data can obscure these corrections, as the drift can 
be corrected with a slightly erroneous exponent. Where 
we have used scaling plots, we do not try to collapse all 
of the data onto a single curve, but keep in mind that 
the finite size corrections give a consistent drift with sys- 
tem size and we therefore tried to optimize the fit to the 
largest systems and near he- 



III. SUMMARY OF RESULTS 

As we are interested in the behavior of the RFIM in the 
thermodynamic limit, we have studied the approach to 
the infinite- volume limit using finite-size scaling analysis 
techniques similar to t VLC jS O || anpliod t o spin glasses and 
other random systems.l^Tl^'1rTr^lr4r'i However, in con- 
trast with ground state studies of spin glasses and other 
random models for which only a single thermodynamic 
phase exists, the results presented riiere give insight into 
the transition between two phases.cJ 



A. Stiffness 

The fundamental difference between an ordered phase 
and a disordered one is the stiffness (or rigidity) of the 
former: the free energy cost of changing one part of a sys- 
tem with respect to another part far away. At a macro- 
scopic level, this free energy cost must be at least of order 
ksT and is usually much larger, diverging as a power of 
the system size. For an Ising ferromagnet, this stiffness is 
provided by the free energy cost of a domain wall which 
scales as its surface area. Thus a natural quantity to 
study for the ground states of the RFIM is the domain 
wall energy. This can be obtained from the difference in 
energy between antiparallel and parallel boundary con- 
ditions imposed on opposite sides of a system of cross 
sectional area L^. A particular combination of these we 
call the stiffness^ which we denote by E. Because of the 
randomness, this energy will be sample dependent and 
there is information to be gleaned from its distribution 
as well as its mean. 

The scaling theory of the putative critical point of the 
random field Ising model predicts that the distribution 
of the stiffness will have a scaling form near the critical 



point: 
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with 6 and v universal exponents, P a universal scaling 
function (which does, however, depend on the shape of 
the sample), and C and K non- universal coefficients. In 
the ferromagnetic phase, the distribution of the stiffness 
will be sharply peaked at long length scales about a mean 
value which grows as a{h)L'^ with cr(/i) the interfacial ten- 
sion. This interfacial tension vanishes as h y' he. In the 
disordered phase, E will typically fall off exponentially 
for system thicknesses, L, much larger than the correla- 
tion length ^{h) ^ {h — hc)~'^ ■ This exponential decay 
of the stiffness with L is confirmed for all values h > he 
examined in our numerical results. 

At the critical point, the distribution of E will be broad 
with both mean and width of order L^ . The exponent 9 
thus characterizes the scaling of the stiffness at the criti- 
cal point. As long as 6 is positive, the basic features of the 
zero-temperature critical point will be stable to thermal 
ffuctuations and the finite-teniperature transition will be 
in the same universality classB 

Our studies of the stiffness are based on computing 
energies for samples periodic in two directions and having 
fixed uniform boundary spins on the other two faces. Our 
results are very consistent with the scaling predictions for 
the larger system sizes, up to 128^, close to the critical 
point, which occurs at 



2.270 ±0.004. 



(3) 



This location for the critical point is consistent with those 
obtained from scaling analyses of the domain wall dimen- 
sion and the magnetization. It is somewhat lower tt 
some previously reported estimates such as he ~ 2.33,1 
but it is consistent with the values he — 2.29 ± 0.04 re- 
ported by Hartmann and NowakiJ, he = 2.26 ± 0.01 re- 
ported by d'Auriac and Sourlas,El and he — 2.28 ± 0.01 
reported by Hartmann and Young.Ea Taking he — 2.270, 
the exponents that give a good scaling fit for the stiffness 
are found to be 



and 



61 ~ 1.49 ±0.03 



i/~ 1.37 ±0.09. 



(4) 



(5) 



The value for 9 is consistent with exact bounds as well 
as with values derived from finite temperature simula- 
tions by applying exponent relations to measured critical 
behavior £3 Note that if the transition had been first or- 
der, one would have expected to find 9 ~ d — 1 = 2, 
with a double peaked distribution of E corresponding 
to "ordered" and "disordered" samples, and an effective 
1/ = 2/d = 2/3; the results we find are far from these. 

From the modified hyperscaling law appropriate i-te. 
transitions governed by zero temperature fixed points pEI 



{d-9)v^2-a 



(6) 



with d the dimension, here equal to three, we predict 
that the specific heat exponent for the finite-temperature 
transition (and for the second derivative of the energy 
with respect to h at zero-temperature near the transition) 
is 



2-(3-6l)i/~-0.07±0.17. 



(7) 



We also fit the sample averaged bond part of the en- 
ergy, Ej, at he to the form Ej ~ ci — C2L'^"'~^'> I ^ to 
more directly obtain a, inspired by the recent approach 
of Hartmann and Youngpj who examined the scaling of 
the derivative dEj /dh. We find a consistent value for a 
using similar methods, although our value disagrees sub- 
stantially with that of Hartmann and Young. We also 
use an extrapolation of Ej^ based on the dimension of 
the domain wall surfaces, which define Ej, to find 

a = -0.01 ±0.09. (8) 

Both values are consistent with experiments,E3 which 
yield a small value of a. 



B. Domain walls 

In addition to the stiffness measurements, we have in- 
vestigated the properties of the domain walls that are 
forced by appropriate changes of the boundary condi- 
tions. In the ferromagnetic phase, we expect that these 
will be flat on large length scales and have area propor- 
tional to L^ . These walls will be rough with a transverse 
width W^ on a scale i. described by the roughness expo- 
nent C,, W ^ i''. But at the critical point, we expect 
the walls will become fractal. The definition of a domain 
wall in the RFIM has ambiguities because some isolated 
clusters of spins — in particular those with anomalously 
strong random fields — are frozen, i.e., unaffected by 
changes in boundary conditions. The identification of 
the bonds that define the domain wall is therefore un- 
certain up to these fixed spins. We use three methods 
for calculating the fractal dimension of the domain walls 
introduced by changes in the boundary conditions; each 
definition has a distinct physical import in a scaling pic- 
ture of the RFIM. 

One method is to determine the surface area of the 
set of spins connected to one face of the sample that 
are unchanged when the spins on the opposite face are 
reversed. This yields a spanning surface of a dimension 
that we denote ds- The second method is based on a 
box counting approach that counts which volumes in a 
system with antiparallel boundary conditions differ from 
both the "up" and the "down" configurations obtained 
from parallel spin boundary conditions; this we denote 
dj, to indicate its role as a measure of the volume locally 
incongruent with these two configurations. 

A third method does not measure a dimension directly 
but rather an energy: the contribution of the exchange 
interactions to the stiffness at the critical point. As will 



be explained later, this fractal dimension, dj, is not ex- 
pected to be an independent exponent. Rather, it is re- 
lated to the others by 



dj = 



1/iy, 



(9) 



a relation obtained by considering the derivative of the 
stiffness with respect to h. It can be seen that our results 
for the exponents are entirely consistent with this scaling 
law. 

The three exponents associated with the fractal dimen- 
sion of the critical domain walls are similar, but perhaps 
not all mutually consistent, given the estimated error 
bars: 



ds = 2.30 ±0.04 
di = 2.24 ±0.03 
dj = 2.18 ±0.03 



(10) 

(11) 
(12) 



As we will discuss, we believe that at least two of these di- 
mensions indeed measure slightly distinct quantities, due 
to frozen spin clusters that are relatively independent of 
boundary conditions. The simplest plausible conjecture 
is that dj = di < ds, though it may be that dj and dj 
are distinct. 



C. Magnetization 

As mentioned above, some previous studies have found 
that the magnetization appears to be discontinuous at 
the transition. Indeed, our data for the magnetization 
as a function of size for various types of boundary con- 
ditions that were chosen so that some favor a ferromag- 
netic state while others favor a disordered state, appear 
to be consistent with the coexistence of three states at. 
the critical point as was found in other recent work.ES 
But, as discussed below, we believe that this conclusion 
is influenced by the nearness of fi/v to zero. Based upon 
the scaling picture and numerical evidence, we will ar- 
gue that, at the critical point there is only a single state 
and the apparent "up" and "down" configurations do not 
correspond to distinct states in the infinite volume limit. 

Using the magnetization data and the best fit critical 
point found from our studies of the stiffness, we can at- 
tempt to extract an estimate for the scaling of the mag- 
netization with system size at the critical point. This 
yields 



f3 



- =0.012 ±0.004 



[magnetization], (13) 



which is inconsistent with zero at the level of three stan- 
dard deviations. The primary uncertainty in our esti- 
mate of P/h' arises from the uncertainty in the value of 
he, as the statistical errors in the sample average |r?T.| are 
relatively small at fixed h. This exponent describes the 
magnetization very well for systems of size 32 < L < 256, 
for a range of he, 2.265 < he < 2.275. 



D. Spin clusters and walls 



The value of 



In spite of the smallness of /3, useful information on the 
decay of spin correlations at the critical point can be ob- 
tained indirectly by studying the statistical properties of 
the domain walls separating connected clusters of paral- 
lel spins. In the ferromagnetic phase, we expect that the 
probability of finding a region of diameter £ that is not 
affected by the boundary conditions decays exponentially 
for £ » e as exp(-C(£/e)''-2). 

Since the system appe ars to be ferromagnetic at the 
critical point, due to |m| being nearly unity for the sys- 
tem sizes studied, we also study clusters of the minority 
spins at /i « he- The clusters are defined hierarchically 
starting from the largest connected cluster of connected 
spins, with the surface of each cluster given by its outer- 
most surface, that is, the set of bonds connecting it to the 
surrounding cluster. The volume of each cluster includes 
that of the fully enclosed subclusters of the opposite sign 
(and their subclusters, if any, etc). But the outer surface 
of a cluster does not include the surfaces of its fully en- 
closed subclusters, whose number scales with the volume 
of the cluster. 

The outer cluster surfaces are found to be fractal with 
mean area a (averaged over clusters and samples) scaling 
with enclosed volume v as 



a{v) 



.fi/d 



with the exponent 



d 



0.755 ±0.008, 



(14) 



(15) 



suggesting a surface fractal dimension d^ ~ 2.26 ± 0.02, 
a value consistent with the domain wall dimension ds . 

Perhaps more interesting is the distribution of the 
number density of spin clusters as a function of their size, 
in particular the probability p{v) that a given site is in a 
minority spin cluster of size of order v — more precisely, 

V 

p{v) = -— Probfsite G cluster of size in (v, v + Sv)]. (16) 
ov 

We find that over the range of sizes studied p{v) ap- 
pears to converge to a small constant value, poo, for 
1 ^ w ^ i^, with periodic boundary conditions. This 
implies that in the limit of an extremely large system, 
any given spin will definitely be in such a "minority spin" 
cluster; indeed, it will typically be within one such large 
cluster which itself will be within a cluster of typical size 
^ p^ larger which itself will be in an even larger cluster, 
etc. This is exactly the type of behavior that gives rise to 
power law decay of spin correlations at a critical point on 
sufficiently long scales, as is explained in Section IX. It 
is consistent with expectations from other observations 
we have made, in particular that the probability that the 
stiffness of a finite sample is exactly zero tends to a non- 
zero constant for large system sizes at the critical point. 



Poo =i 0.0019 ± 0.0004 (17) 

that we findEEl yields an estimate for 

/?/i/ = 2dpoo =i 0.011 ±0.003 [cluster] (18) 



consistent with Eq. (|13|). This exponent controls the de- 
cay of the typical magnetization with system size at the 
critical point: 



i{K) ^ L-^'\ 



(19) 



For L = 128, this only gives a reduction factor of 0.94 
from the magnetization of a small system and is consis- 
tent with our magnetization data. Note that with this 
estimate, one would need to go to system sizes of order 

10^^ ~ L ~ lO'^^ to see a factor of two reduction in the 
magnetization at the critical point! 

For the magnetization in the ferromagnetic phase, us- 
ing this calculation of P/v and consistent with the finite 
size scaling of the magnetization, we expect conventional 
behavior with 



but with 



mr^ {he- hf 



/?- 0.017 ±0.005 



(20) 



(21) 



— far smaller than for any other known system with 
the exception of the one-dimensional Ising model with 
long-range l/r^ interactions which has a critical transit, 
tion with a discontinuous magnetization, i.e. /3 = 0.E2I 
This small value is near the value /3 — 0.02 suggested 
from numerical rerprmalization group calculations on a 
hierarchical lattice.E2l The numerical value we find is con- 
sistent with several previous studies: for example, Hart- 
mann and Nqssiak determine^ = 0.02 ±0.01, using exact 
ground^tatesEl, Swift, et al^ find (3 = 0.025±0.015, and 
RiegeiO found /3 « at positive temperature, but with- 
out any latent heat or multipeak structure in the mag- 
netization distribution, suggesting a second order transi- 
tion. However, we can more clearly exclude [3 /v = as 
a possibility by making use of connections between the 
value of (3/v and the statistics of spin clusters. 

One question that naturally arises concerns the struc- 
ture of spin clusters for small |m| near the critical point, 
where the sample no longer appears ferromagnetic. Is 
there a possibility of percolation of both up and down 

spins, when \h — hc\ ^ 10~(^^='=^) in large enough sam- 
ples? For fixed ± or — boundaries, as /i ^ he, \m\ be- 
comes small, but the minority and majority spins are not 
independent. Hence, even though the density of + and — 
spins becomes almost equal, the minority spins are large 
clusters embedded within the matrix of majority spins, 
so that only one sign of spin percolates in the disordered 
phase even close enough to the critical point that the 
magnetization is very small and the density of minority 



spins almost one half. Exactly at the critical point in an 
infinite sample, or where ^ > L in a large finite sample, 
the long length scale characterization of the spin con- 
figuration will be rather different than in either phase; 
these differences motivated some aspects of the present 
numerical study. 



E. Number of states 

To study the RFIM phases and transition in more de- 
tail, we have analyzed the inffuence of boundary condi- 
tions on a window of size w in the center of a sample as 
the sample, size L diverges. As made clear by Newman 
and SteinEil, the character of the thermodynamic limit of 
the ground states can be investigated by studying such 
windows. Our numerical computations strongly support 
the picture of a small number of ground states — two 
in the ordered phase and one in the disord er e d .phase — 
consistent with the simple scaling scenario .c1Ie3 

Nevertheless, because f3/iy is so small, at the critical 
point it it is difficult to use numerics to directly distin- 
guish between two scenarios: (A) two coexisting states, 
as in the ferromagnetic phase; or (B) a single state, 
with interior spins unaffected by boundary conditions as 
i — > oo. up were exactly zero, as in (A), then the proba- 
bility q that boundary conditions could affect spins in the 
center in ways other than the apparent "up" and "dos 
phases would decay as a power of the system size^ 
q r^ L'^='~'^, where dg is the fractal dimension of domain 
walls. In scenario (B), a similar power law scaling is ex- 
pected, with q ^ i^ds-p/v-d^ where the change in expo- 
nent reflects the freezing of spins or, equivalently, decay of 
magnetization, as i — > cxd. As will be argued below, the 
simplest expectation is that the exponent dj — d^ — P/v, 
so that q ^ L"^'^'^. This is consistent with the assump- 
tion of only one state at criticality. 



IV. OUTLINE 

The remainder of the paper gives the details of the nu- 
merical results and related scaling arguments. Table 
is a summary of the numerical values of the exponents. 
In Sec. M, we describe how the stiffness is computed and 
demonstrate that its scaling is quite consistent with a 
"conventional" second order phase transition. To aid in 
developing understanding, we study how the probability 
that the stiffness is exactly zero depends on the sample 
shape. Sec. VI presents the three methods that we em- 



ploy to compute the dimension of the domain walls gen- 
erated by comparing different boundary conditions (the 
same comparison used when calculating the stiffness.) 
The methods differ somewhat in how they count regions 
of "frozen" spins that are not affected by boundary con- 
ditions. In the subsequent section (Sec. VII), we report 



the scaling of these distributions are quite consistent with 
a single value of he (and also consistent with the methods 
of finding he in other sections.) Our study of the scaling 
of the surfaces o f spin clusters with their volumes is sum- 
marized in Sec. VIII, Besides giving a fractal dimension 
d1 consistent with the domain wall dimension dg, these 
computations can be used to separately infer 13/ v, given 
an understanding of magnetization and correlation func- 
tions based upon a domain wall picture. Our estimates 
for the singular behavior of the specific heat are included 
in this section. The general scaling picture that connects 
these results is reviewed in more detail in Sec. IX. In Sec. 



M, we report results of how the spin configurations de- 
pend on sample size and boundary conditions for a fixed 
disorder realization. These results are consistent with a 
single transition separating a (large /i, disordered) phase 
with a single thermodynamic limit from a (small /i, or- 
dered) phase with two distinct thermodynamic limits. In 
the Summary (Sec. XI), we review the scenario for the 
transition that is consistent with the numerical results 
and contrast this scenario with alternate pictures. 



V. STIFFNESS AND SCALING 

To establish the location and nature of the transition, 
we first focus on the stiffness of the system. In an or- 
dered Ising phase, the (free) energy of a domain wall 
across a system of size L'^ will be E « aL'^~^ with a the 
interfacial tension. At an ordinary first order transition, 
the interfacial tension is discontinuous at the transition, 
while near a second order transition, it goes smoothly to 
zero. For a zero-temperature transition, the interfacial 
tension vanishes with a variant of Widom scalingj 



-h\(d-l- 



(hc - h) 



(22) 



results on the magnetization m near he- Though the dis- 
tribution of m depends strongly on boundary conditions. 



To probe the stiffness of a random system takes some 
care. For a random field system, there is no exact sym- 
metry between the up and the down spins but only a 
statistical symmetry of the distribution of the random 
fields. Thus, for example, for a given sample in the or- 
dered phase, the energy of the up state — obtained with 
up boundary conditions (BCs) — will differ from that of 
the down state (obtained from down BCs), by a random 
amount of order vV arising from the differing effects of 
the random fields on the two states. In order to compute 
the interfacial energy it is useful to subtract as much as 
possible of this random "bulk" energy so as to be left with 
a quantity that is as-close-as-possible to an "interfacial" 
energy. 



A. Definition of the stiffness 

To obtain the stiffness of a sample, we compute the 
symmetrized energy difference between antiparallel and 
parallel boundary conditions. This is computed from the 



TABLE I: Table of numerical estimates. The exponent or constant name, computed value, primary method for inferring the 
value, section discussed, and most relevant figure are listed. 



Symbol Value 



Definition and data used 



e 1.49 ±0.03 



Scaling of stiffness at he, violation of hyperscaling. 

Found from scaling of stiffness with L and h — he, see Sec. [V| and Fig. 



he 2.270 ± 0.004 Critical value of the random field. 

Determined from constancy in Pq, probability of zero stiffness (see Fig. ^ and Sec. |v|), 
and consistent with estimates from convergence of effective dimension estimates ds,i,j, 
scaling of A„2 peak locations with L, scaling of \m\ with L (Fig. M), 
and window change probabilities (Fig. p9[). 



v 1.37 ±0.09 



Correlation length exponent. 

Found from scaling of the stiffness with L Sec. M and Fig. 

Consistent with scaling of A„2 peak locations with L. 



with he fixed by Pq measurements. 



C 0.66 ±0.03 



Roughness of domain walls in the ferromagnetic phase. 

Found using anisotropic scaling (Fig. nsl) and effective exponent in L^' "^ x L^ samples. 

See Sec. Ivf^. 



ds 2.30 ± 0.04 



Fractal dimension of connected domain wall at h — he- 
Found from surface of W++,^ , as shown in Fig. [ll[ 

See Fig, ^and Sec. [VIB|. 



di 2.24 ± 0.03 



"Incongruent" fractal dimension of domain wall at criticality. 

Box counting of incongruent volumes (disconnected wall). See Fig. 15 

Consistent with scaling of state overlap probabilities shown in Fig. 29 



and Sec. VfD 



dj 2.18 ±0.03 



Energy "fractal dimension" at h = he. 

Found from the exchange part, Ej, of the stiffness. 

See Fig. pi|and Sec. IviE . 



d": 2.27 ±0.02 



Fractal dimension of the surface of spin clusters. 
See Fig. p2|and Sec. |VIIIA , 



Poc 0.0019 ± 0.0004 Probability per scale e of crossing a spin cluster surface a.t h = he 
See Sec. |VmB| and Figs. p3|, psf and |2^. 



iS/f 0.011 ± 0.003 Ratio of magnetization exponent to i/. 

Determined from poo and consistent with scaling of \m\ vs. L at criticality. 
See Fig, ^and Sec. ItXAJ 



13 0.017 ± 0.005 Magnetization exponent, found from /S/f and i/. 



{a-l)/iy -0.74 ±0.02 



Combination of heat capacity exponent a and v 
Found using value for dg and Eg. (p6|). 



a -0.01 ± 0.09 



Heat capacity exponent, found using Eq. (p6|) and u. 
Consistent with modified hyperscaling Eq. (Bl) and the value a - 
from a fit to the bond energy density Ej{L) at he and v. 



'0.12 ±0.12 found 



ground state energies for four different boundary con- 
ditions on a given sample, denoted +±, H — , — h and 

. These correspond to fixing the spins to have values 

s = +lors = — Ion the left or right sides while imposing 
periodic boundary conditions in the other two directions. 
For example, H — has spins fixed to +1 on the left and 
tO|-^l on the right. The interface energy is then defined 



(i?+_ + E^, 



E, 



E-)/2. 



(23) 



Note that the average over samples of S will be the same 
as that of E^w = E^ — -£'++■ Studying E, however, re- 
duces the effects of energy changes near the boundaries 



that are caused by the differing boundary conditions: in 
E, each boundary condition on each side appears twice 
but with opposite signs so that these effects cancel. This 
cancellation will be most pronounced well into the disor- 
dered phase. 

In the disordered phase, the boundary conditions typ- 
ically only affect layers near the boundaries with thick- 
ness of order the correlation length ^; deep in the inte- 
rior (for system sizes L ^ ^) the spins will be frozen, 
completely unaffected by the boundary conditions. The 
average energy of the boundary layers will, because of 
the statistical symmetry, be independent of whether the 
boundary conditions are plus or minus. But there will 



be a random part of the boundary energy, with magni- 
tude of order \/ L'^~^^ that is sensitive to the boundary 
condition. Thus in three dimensions, i^+H^ will typically 
be of order L even in the disordered phase. In contrast, 
the stiffness E will typically be exactly zero because of 
the cancellation of the boundary energies and the con- 
comitant frozen interior which blocks any knowledge of 
the spins near one face about the boundary conditions 
on the opposite face. In general, the distribution for S 
contains a (5-function contribution with some weight 



(a) 



Po = Prob[E = 0]. 



(24) 



Sample configurations from simulations are illustrated 
in Fig. gj with two-dimensional slices shown. Part (a) 
of the figure illustrates a situation somewhat into the 
disordered phase in which the left and right boundaries 
are effectively decoupled as discussed above. The frozen 
spins, those that are the same with all four boundary 
conditions, are indicated by dark or white squares in part 
(c) of the figure, while those that are affected by the BCs, 
the controllable spins, are indicated by gray squares. 

The behavior in the ordered phase is quite different as 
can be seen in the parts (b) and (d) of Fig. ||. In this 

case, the difference between the , the H — and the 

H — h boundary conditions can be well characterized by a 
+ |— domain wall that has a minimum energy position 
somewhat to the left of the center. Similarly the differ- 
ence between the , the — \- and the -I — h boundary 

conditions is characterized by a — 1+ domain wall whose 
minimum energy position is somewhat to the right of the 
center. The stiffness of this sample will thus be half the 
sum of energies of the two types of walls plus the energy 
of the random fields (here predominantly negative) in the 
region between the two favored positions of the walls; the 
contribution of the random fields in this region will not 
cancel. 

This picture yields a stiffness in the ordered phase with 
a mean of order L^^^ — L^ and variations around this 
mean of order L'^''^ — L^ , the variations being dominated 
by the random fields in between the positions of the two 
types of walls. In the ordered phase for h < he, Po ^ — 
apparently exponentially fast or faster in L — as L — *■ oo. 

At the critical point, the behavior is qualitatively like 
that in the ordered phase. But here the energy cost of 
the interface is much lower, the interface itself is fractal, 
and, in the regions of controllable spins that are otherwise 
flipped by the changing boundary conditions, there are 
large frozen unflipped "holes" . 



B. Statistics of the stiffness 

The stiffness S was determined by finding the ground 
states for a single sample subject to each of the four 
boundary conditions -|--|-, , — h and H — . By study- 
ing many samples, we computed the distribution of E for 
various system sizes and random field strengths. In par- 
ticular, we computed both the probability Pq that E = 





L = 40, h = 2.800 



(b) 





L = 40, h = 2.200 



h^. 



L = 40, h = 2.800 




L = 40, h = 2.200 



FIG. 1; Pictures of planar slices (z = ) of configurations, 
for fields (a) h = 2.8 and (b) h — 2.2, in samples of size 
40^. The slices shown at each h are the four combinations 

of boundary conditions (top left), — h (top right), H — 

(bottom left), and -I— I- (bottom right), where the left and 
right faces (in the x direction) are fixed -I- or — and periodic 
boundary conditions are in effect for the y (up/down) and 
z (out of the page) directions. The dark squares indicate 

an up spin at that location in the slice. The and ++ 

visualizations for ft = 2.2 show the presence of minority spin 
"bubbles" embedded in the bulk. A summary of the effect of 
the boundary conditions for h = 2.8 and h — 2.2 are shown 
in parts (c) and (d), respectively. Dark and light squares 
correspond to up and down spins, respectively, that are frozen, 
i.e., invariant under this set of boundary conditions. The 
gray controllable spins can be modified by choosing among 
the four boundary conditions. For h = 2.8, the gray volume 
is composed of two unconnected regions anchored on the two 
controlled boundaries, so that the stiffness E = 0. In contrast, 
a,t h — 2.2, in the sample shown, the gray region connects the 
two sides and E 7^ 0. 



as well as the mean stiffness E, denoting, as usual, av- 
erages over the randomness by overlines. For the bulk 
of the computations, ground states were found for cubic 
samples of size L^ and anisotropic samples of size 2L x L^ 
with the length along a (100) axis of the lattice perpen- 
dicular to the controlled faces, the a;-direction, being 2L. 
To check that the results were not artificially influenced 
by lattice orientation effects, we also computed values of 
S for two types of samples whose controlled faces are 
L X L rhombi, with L and 4L layers, respectively, sepa- 
rating the two faces along the (111) direction. Note that 
such (111) layers are separated by a distance of 1/V3, 
rather than the distance of 1 that separates the (100) 
layers. As the lattice in both of these geometries is cu- 
bic and the ferromagnetic couplings are the same, the he 
found should be the same in the two orientations. 

If the transition is second order, the mean interface 
energy should scale as 



T" 
A 



i:KCL'^S[L^^''{h~hc)K' 



(25) 



where the exponent 9 sets the scaling of the energy at 
the critical point, f is the correlation length exponent, S 
a universal scaling function that depends on the shape of 
the sample, and K and C are non-universal cocfhcients. 
Using this scaling form and varying 6, v and h^ yields a 
good collapse of the data, as shown in Fig. Hfor the (100) 
samples. By varying the exponents in the scaling plot, we 
estimate the values h^ = 2.27 ±0.01, 9 = 1.50 ±0.08 and 
V = 1.35±0.20. Using the data for Pq to fix h^ = 2.270 as 
discussed below, gives 9 = 1.49±0.02 and v = 1.37±0.09. 
The excellent collapse of the data strongly support the 
conclusion that the phase transition is second order. 

The value of 9 is in quantitative agreement with resuks. 
from Monte Carlo simulations at finite temperature,t3 
which found 9 = 1.53 ±0.10. It is also within the bounds 
determined from various arguments: 



2 i/ - - 2' 



(26) 



the lower bound-jarising from scaling laws and a rigor- 
ous inequalityao. The upper bound follows from the 
observation that any larger value of 9 would imply that 
the system would be stable — by the argument of Imry 
and MatJ — to an increase of the random field and thus 
should not be at the critical point. Since (3/v is extremely 
small, we expect that the true value of 9 should be just 
slightly below 3/2. This is to be contrasted with the 
"dimensional reduction" result predicted to obtain to all 
orders in a d = 6 — e expansion oi 9 — 2 (but see recent 
claims in Ref. pSJ ). 

The correlation expoiaent v must be no smaller than 
2/d in random systems.ES Our result easily satisfies this 
bound. Indeed, it is substantially larger than this lower 
bound and even more so than the mean-field value of one 
half; this is presumably associated with proximity to the 
lower critical dimension of d^ = 2. 
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FIG. 2: Scaling plot for the sample averaged stiffness E for 
(a) isotropic samples of volume L^ and (b) anisotropic samples 
of volume 2Lx L^, with the longest axis being the direction in 
which the boundary conditions are varied. Note that the ver- 
tical scales differ. The stiffness is calculated by the symmetric 
comparison of four ground state energies: the energies for the 
four choices of spin up and spin down boundary conditions 
on the left and right sides and with periodic boundary condi- 
tions in the other two directions. The fit shown is for energy 
exponent 6 = 1.49, correlation length exponent v — 1.37 and 
critical value he = 2.270. This scaling is consistent within 
errors, except for the L = 16 isotropic samples. Statistical 
(la) error bars are shown. 



1. Stiffness in the disordered phase 

Fig. S shows the dependence of the mean stiffness on 
the linear dimension L for the L x L^ samples. The de- 
cay of the stiffness is well fit by a decaying exponential 
S ~ exp(— L/^s), for h > he and L » ^s (roughly 
when E < 0.2.) The correlation length can be inferred 
from the fits. The values for ^s obtained from the 2L x L? 
samples, using a similar plot, are in agreement with those 
from the cubic sample to within 10% for each h. The val- 
ues of the correlation lengths ^s found are consistent with 
a divergence of 5e ~ (/i — ft.c)^^'^^^'^, taking h^ = 2.27, 
consistent with our other determinations of v, though the 
data is not very near he- It may be possible to make a 
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FIG. 3: Plot of the decay of the mean stiffness E with L, 
in the disordered phase. The hnes are fits to the data for 
S < 0.2 of the form E ~ exp(— L/^e). In conjunction with 
similar fits for 2L x L samples, which allow us to estimate 
errors from finite size fitting, we find the values ^e = 26 ± 4, 
7 ± 1 and 4.0 ± 0.4 for h = 2.4, 2.6, and 2.8, respectively. 



more accurate determination of ly by more careful calcu- 
lations using h values somewhat nearer to he- 

For the 2L x L^ samples, Pq, the probability of the 
stiffness being zero, can be appreciable for accessible L 
and h near h^- The data for fixed h > he are consistent 
with Pq approaching one exponentially with L, although 
other forms cannot be ruled out. 



2. Stiffness at criticality 



At the critical point, a non-trivial scaling function 



Pc(S/CL'') = P(S/CL^O) 



(27) 



for the distribution of E, would suggest that Pq should 
approach a finite fixed point value. In Fig. H(a) a plot 
of Po is shown as a function of L for various h. Observe 
that opposite boundaries are almost always coupled at 
he in the cubic samples: 



^cubi( 



{he) -0.04 ±0.01 



(28) 



This value for Pq is so small that to verify that Pq indeed 
approaches a non-zero constant at the transition, we also 
performed simulations for anisotropic samples of various 
shapes. In general, we expect that the distribution of S 
at the critical point, Pc^ will depend on the shape of the 
sample with long thin samples typically yielding lower 
stiffness and a higher probability of the stiffness vanishing 
than short fat ones of the same cross section. 

For rectilinear samples of dimensions 2L x L^, with the 
controlled boundaries at opposite ends of the long (100) 
axis, we find that Pq approaches a value well away from 
zero, Po = 0.298 ± 0.005, at /i = 2.270 (Fig. |(b)). For 
rhomboidal samples with L^ spins consisting of L lay- 
ers and a length along the (111) control axis of P/a/3, 



(a) LxL 



n 



• h = 2.000 
Ah = 2.240 
<h = 2.270 
■ h = 2.300 
Oh = 2.400 

♦ h = 2.600 



0.05- 



n 



I 



-n 



0.00 L 
0.50 r 

0.40 



10 



20 



50 



100 



'TT 



_o 



(b) 2L X L 







• 


• 


^ 






0.30 


4ft 


^ 








- 


^ 


-<l- 






Oh =1.800 






^ 


V 




• h = 2.000 







"O 




A: 


V 


D h = 2.200 




0.20 


- 


o 




A; 


Ah = 2.240 

Vh = 2.255 


- 












-^ h = 2.270 






# 








Oh = 2.285 




0.10 








o 


■ h = 2.300 


- 




-o 


^ 






Oh = 2.400 






o 


* 


A 


, 1 


0.00 jQ 


20 


L 


50 * 


100 


50 


1 




1 1— 


— \ — ' — ' — ' 


1 1 



0.40- 



0.30- 



0.20^ 



0.10- 



0.00 L 



1/2 1/2 2 

(c)(lll),4L/3 x3 r 



Ah = 2.240 
-4 h = 2.270 
■ h = 2.300 



10 



20 



50 



100 



FIG. 4: Plot of the probability Po{h,L) that the stiff- 
ness E is zero, for (a) isotropic samples of volume L^ and 
(b) anisotropic samples of volume 2L x L^, with the longest 
axis being the direction in which the boundary conditions are 
varied, and (c) anisotropic samples of volume 4L , with the 
boundary faces in the (111) plane. For all sample shapes, the 
convergence to a fixed value of Po as L ^ oo for h — 2.270 
suggests the location of the critical point. The solid lines con- 
nect the points for h = 2.270 to demonstrate convergence of 
Po to a constant, within statistical errors. As Po is very nearly 
zero for isotropic samples (Po(2.27, oo) ~ 0.04, if the appar- 
ent convergence holds at large L), the errors in determining he 
are larger. From the 2L x L^ anisotropic samples, where the 
apparent extrapolation is Po(2.27, oo) = 0.298 ± 0.05, he = 
2.270 ± 0.004. For the (111) oriented samples, with volume 
of 41// v^ X y/SL^ (layer separation x layer area), the data is 
also consistent with he = 2.270, with Po = 0.23 ± 0.01. Sep- 
arate results, not shown, for (100) samples of shape 4L x L^ 
give a value of Po(2.27, 16 < L < 64) = 0.79 ± 0.02. 
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FIG. 5: Probability density p{a) of the values of the inter- 
facial energy density a = L~^Y^ in the ordered phase with 
h = 1.60 < he- As the sample size grows larger, the relative 
sample-to-sample variations of a decrease, consistent with an 
approach to a 5- function at the mean value a{h = 1.6) « 0.69, 
for both the L^ and 2L x L^ samples. 



FIG. 6: Probability density p(S) of the scaled non-zero val- 
ues of EL"* with e = 1.49, for h = 2.270 ^ K. A 5-function 
at E = with weight Pq ^ 0.04 [^ 0.298) for the L^ (2L x L'^) 
samples is not shown. The distributions for samples with lin- 
ear dimensions L greater than 16 are statistically consistent 
with a fixed point distribution for E, with a characteristic 



scale Eo 
samples. 



and a form dependent on the shape of the 



Pq is not distinguishable from zero. However, for longer 
rhoniboidal samples with AL'^ spins consisting of 4L lay- 
ers and a length AL/y/S along the control axis, we find 
Po = 0.21±0.01at/i = 2.270forl6< L< 64 (Fig.|(c)). 
Imposing convergence of Pq to a fixed (non-trivial) value 
as L -^ oo gives a critical value of he — 2.270. These 
anisotropic rectilinear and long rhomboidal samples yield 
our most precise estimate for he, Eq.(||). 

We should expect Pq to be a smooth function of the 
shape; the fact that it is far from zero in samples with 
aspect ratio of order two lends strong support to the con- 
jecture that it will be non-zero for any shape. The ob- 
servation that it is small for cubical samples is related, 
as will be explained below, to the smallness of /3. 



3. Comparison of distributions for E 



The complete probability distributions for S at vari- 
ous values of h are plotted in Figs. || through Fig. || for 
both the cubical and the elongated (100) samples. For 
h = 1.6, the distribution of E appears to approach a nar- 
row distribution about S « (1.39)L^, regardless of the 
sample shape; this is as expected for the ordered phase. 
For h = 2.27, the critical point, the distribution obeys the 
simple scaling form of Eq. (||) for both the isotropic and 
anisotropic samples but with a different scaling functions 
for each of the two shapes. Fig. g shows the integrated 
probability distributions for /i = 2.40. As L increases, 
the mean interfacial energy decreases approximately ex- 
ponentially, and Pq approaches 1. 
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FIG 



7: Scaling plot for the cumulative distribution /(E) = 
J^ dT,'p{T,') of the stiffness, for h = 2.270 ^ he. The stiffness 
E has been scaled by the energy scale L", with 8 = 1.49. The 
labels indicate the sample shapes (2L x L^ and L x L^) for 
each set of curves. For each sample shape, four sample sizes 
are plotted (16 x 8^ ^ 128 x 64^ and 16^ -^ 128^.) At the 
resolution shown, the scaled curves are nearly independent of 
L. The intercept at E = corresponds to 7(0) = Po, the 
probability of a sample having zero stiffness. As in Fig. h, the 
curves converge to a fixed point distribution. 



VI. GEOMETRY OF DOMAIN WALLS 

In addition to the scaling properties of the energies of 
domain walls, we are also interested in their geometrical 
properties. These properties are expected to be related to 
the properties of the surfaces of spin clusters that are ei- 
ther frozen or induced by bulk perturbations (as opposed 
to boundary perturbations), to the effects of boundary 
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FIG. 8: Cumulative distribution /(ct) = J da' p{a') for /i = 
2.400 > he- As the sample size increases, Po, given by the 
intercept of the curves at E = 0, increases, and the typical 
non-zero values of cr = EL~ rapidly decrease. From Fig. pi 
the length scale for the decay of the stiffness is ^e = 26 ± 4. 
This length is comparable to the mid-range system sizes here. 
Note that Pq rises more quickly and the typical non-zero o 
decays more rapidly with L in the anisotropic samples. 



conditions on the deep interior of a sample, and to the 
general scaling picture for the transition. 

In the ferromagnetic phase, the interfacial tension a is 
positive. The domain walls will appear flatter and flat- 
ter on large length scales with surface area proportional 
to i^ but nevertheless divergent roughness characterized 
by a roughness exponent C, = 2/3 and raBdam-eaergy 
variations that scale as i^^ with Oi = 4/3.BBBEI 

At the critical point, the interfacial tension of the walls 
vanishes. Thus we should not except them to be flat even 
on large scales; the natural expectation is that they will 
be fractal with surface area scaling as 



A^L" 



(29) 



with ds a fractal dimension. One might expect, a priori^ 
that the exponent d^ would be an independent exponent 
as it is not obviously related to 9, (3, and v. For the 
simple scaling scenario to obtain, we expect ds to be in 
the range: 



A. Frozen spin regions 

To study interfaces we would like to compare the spin 
configurations found using the boundary conditions 4—1-, 

H — , — h, and as discussed in the previous section. 

But in random field systems, there is an intrinsic dif- 
ficulty associated with defining an interface: this arises 
from the presence of frozen regions which are not affected 
by changing from H — \- boundary conditions to bound- 
ary conditions and thus are unaffected by any chanmes. 
in the boundary conditions on the controlled surfaces .clI 
With mixed boundary conditions, say H — , the interface 
between the region that is like the "up" (-I--I-) state and 

the region that is like the "down" ( ) state can pass 

along the boundary of the frozen regions. Are we to count 
such sections as truly part of the interface? Or should we 
exclude the frozen regions from the system and think of 
the interface as bisecting only the remaining controllable 
regions? 

We are thus led to consider several methods for mea- 
suring the surface area of "interfaces" anticipating that 
we might obtain results which depend on the definition. 
For these considerations, it is useful to refer to Fig. ||, 
which is a sketch of what might happen when S = 0, and 
Fig. |lO|, which is a representation of what might happen 
for E ^ 0. 

Configurations are shown for each of the four bound- 
ary combinations: the circles enclose regions of "frozen 
spins" — those that are constant under all four BCs 
— with solid lines indicating broken (unsatisfied) bonds. 
The dashed lines indicate the location of a frozen clus- 
ter embedded in a set of like spins. The interior of the 
configurations in Fig. |l^ are also frozen. Note that the 
frozen spin regions can contain nested subclusters of al- 
ternating spins. In Fig. 0, spins outside of the frozen 
spin regions can be either -t- or — , depending on the EC 
combination. These two figures are caricatures of config- 
urations such as those shown in Fig. |l|. Note that Fig. g 
docs not show all of the possibilities. Also, these pictures 
are two-dimensional slices, which hides the possibilities 
of regions having three dimensional "handles" and mini- 
mizes the potential role of simultaneous percolation of -|- 
and — spins in some regions. 



B. Surface exponent ds 



d-1 < d. < d. 



(30) 



If the transition were first order, one would expect dg — 
d—1 as in the ferromagnetic phase (more precisely, some 
fraction of samples at the transition would show such 
interfaces). If, at the other extreme, it were found that 
ds = d, this would mean that the "walls" would be space- 
filling (up to possible logarithmic factors) ; this would cast 
doubt on the overall scaliug scenario for excitations, etc, 
near the phase transition.0 



The first method of defining an interface uses just two 
different boundary conditions, for example, the -I — to 
++ comparison. This change in boundary conditions 
causes a connected set of spins anchored to the right face 
to flip from up to down along with the forced right bound- 
ary spins when the H — \- boundaries are replaced with H — . 

This set of changing sites, which we denote C^ _++, has 

a bounding surface — indicated by the heavy and light 
lines in Fig. O, for the spin conflgurations of Fig. pi But 
some of this boundary will surround islands of fixed up 
spins (some of which themselves have down spin inclu- 
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FIG. 9: Schematics of the spin configurations for four differ- 
ent boundary condition combinations, for a case with E 7^ 0. 
Here, there is a set of controllable spins, connected across the 
sample, that can be either + or — , depending on the boundary 
conditions. These are the majority of the spins in the figure 
shown. The frozen spins are those that are constant under 

the four boundary conditions , — h, H — , and ++; these 

are indicated here by the circular regions. Solid lines separate 
spins of opposite sign, while the dashed lines indicate frozen 
islands that are of the same sign as the surrounding spins. 
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FIG. 10: Schematics of the spin configurations for four differ- 
ent boundary condition combinations, for a case with E = 0. 
In addition to the frozen islands, shown as circles as in Fig. bl 
there is a set of frozen interior spins that spans the sample in 
the directions perpendicular to the horizontal (control) axis. 
Conventions for solid and dashed lines are as in Fig. BI The 
surfaces used to measure ds are the two surfaces of the frozen 
interior, but the measure used to compute di is zero, as long 
as the boxes have side B smaller than the size of the frozen 
interior. Also zero is the exchange stiffness Ej, as each bond 
that is broken in both the H — and — \- configurations is also 

broken in both the and -I— |- configurations, while bonds 

that are broken exactly once under one of the two antiparal- 
lel BCs is likewise broken exactly once under parallel BCs, so 
that all broken bonds cancel in the signed sum that defines 
Ej. 
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FIG. 11: Schematics of the definitions of domain wall mea- 
sures, based on the configurations of Fig. H. (a) The heavy 
solid lines indicate the boundaries used to define the domain 
walls for the calculation of the fractal dimension da of the 
spanning wall obtained by comparing the -I-— and ++ config- 
urations. The region of changed spins connected to the right 

face is C-i ,++, which has both the heavy and light lines as 

boundary, while the unchanged connected region anchored on 

the left face, with the single solid line as boundary, is W^ ,++. 

(b) Boxes used for determining dj, the dimension of the lo- 
cally incongruent regions. The number of boxes of side B 
in which the H — configuration differs from both the ++ and 

configurations scales as L'^' . The broken bonds around 

the frozen islands in the ++ or configurations are not 

counted, (c) The signed sum of broken bonds that defines 
Ej, the exchange contribution to the stiffness E. Solid lines 
indicate positive contributions and the long dashed lines in- 
dicate negative contributions. 



sions) that are disconnected from both controlled faces. 
The number of such islands (hght circles in Fig. [ij) will 

scale with the volume of the C-| ,++ region and their 

boundaries will contribute an amount of order this vol- 
ume to the surface area of C-| _++. This internal contri- 
bution to the surface area will, on average, be dominant 
in large systems when h < he, but it is clearly not prop- 



erly part of the domain wall. 

What we are interested in is the part of the bound- 
ary of C^ ,++ which interfaces with the other "half of 

the system. One way to define a domain wall is thus 
to start at the unmodified left face and find the set of 
spins connected to this face that do not change when the 
boundary conditions on the right face are changed; this 
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set, which we denote lA^ ,_|-+, is simply connected, i.e. 

it has no interior holes, although it could have handles. 
The surface of the set ZY+_ ^++ is just its interface with the 

set C^ ^++ that flips. This hi — C interface, which spans 

the whole cross-section with no holes and thus includes 
some boundary of frozen regions, is our first definition of 
a domain wall of interest. 

Averaging over samples at fixed h gives a mean surface 
area of this U — C domain wall, A{h, L). (For these and 
related studies, we used 5 x 10'^ to 20 x 10'^ samples for 
smaller sample sizes, 8^ through 64"^, and 300 to 5 x 10'^ 
samples for the largest sample size, 128^.) Estimates of 
the dimension of these surfaces, dg, can be obtained from 
the discrete logarithmic derivative, 

ds(h,L) = ln[A(/i, V2L)/A(/i,i/V2)]/ln(2). (31) 

A plot of these estimates, with statistical errors, are 
shown in Fig. ^ The estimates for the case oi h — 
2.27 ~ he appear to approach a fixed value, d^, as 
L -^ oo, while ds ^ 2 for /i < /ic, as expected. For 
h > he, the apparent exponent either starts at ds > dg 
and falls, or first rises before dropping with L. This be- 
havior presumably arises for L <C C, where the growing 
volume allows for larger surface area, while for L ^ £^, 
the domain walls become confined to a distance less than 
^ from the right and left faces of the sample and thus 
effectively become two dimensional. From this plot and 
the results for the (111) orientation (Fig. QTh, we estimate 



2.30 ±0.04, 



(32) 



where systematic errors due to finite size effects and un- 
certainty in he dominate the statistical uncertainties. 



C. Roughness in the ferromagnetic phase 

We have verified that the surface roughness of the 
non-fractal domain walls in the ferrom a | gae1jic . ph ase are 
consistent with theoretical expectations.EZl'E3oc3 Specif- 
ically, we calculated the "height" of the surfaces — devi- 
ation from flat — in anisotropic samples of shape X x L^, 
with the outer two layers in the x-direction fixed to be 
-|- or — and, as before, the sample periodic in the y and 
z directions. Again, to reduce lattice artifacts, we use 
samples whose "x" faces are oriented in either the (100) 
or (111) directions. 

As overhangs are possible in these interfaces, it is nec- 
essary to define carefully the "height" function, u{y,z),: 
for a given y and z coordinates, we use twice the average 

of the a;-coordinates of the set of spins in U^ .__ ; in 

the absence of overhangs, this gives the desired surface 
height . The sam ple averaged rms width W is defined by 
W^ = [u^] — [u]^, where the square brackets indicate the 
average of u{y, z) over the y—z coordinates of the sample. 
Simple scaling in the ordered phase suggests that 




FIG. 12; Effective dimensions ds(h, L) obtained from a loga- 
rithmic derivative of the surface area with respect to L. These 
are used to estimate the fractal dimension ds = ds{hc,oo). 
The values ds (h, L) are calculated from the surface of the con- 
nected set of spins rooted at one face that is unchanged when 
the spins on the opposite face of the sample are flipped. The 
scaling of the area of this surface with L yields the estimates 
shown, via Eq. (pll). The error bars represent la statistical 
uncertainties. The values converge to ds = 2.30 ± 0.04 for h 
near 2.27 « he, with the error reflecting the uncertainty in he 
and the estimated magnitude of finite size corrections. The 
lines connect data points with the same h. 
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FIG. 13: Scaling plot for the roughness of a forced interface 
in the ferromagnetically ordered phase as a function of the 
aspect ratio of an X x L x L sample. For the values of h 
shown here with h < he, the width of the interface scales as 
W ^ L'' with the best fit (^ = 0.64(3), comparable to the 
expected exact result (^ — 2/3. The statistical la error bars 
are 1/5 of the symbol sizes or less. 



for large values of X and L, with T a geometry dependent 
function. We find that using C, — 0.64 ± 0.03, consistent 
with the expected valueE3 C — 2/3, describes the data 
fairly well, as seen in Fig. |l3|. 

The convergence of the roughness of the interface to 
its asymptotic form is made more apparent by defining 
an effective scale dependent roughness exponent 



W = L'^T{h,X/L'^), 



(33) 



C{h, (LiL2)i/2) = ln{W2/Wi)/ HL2/L1 



(34) 
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FIG. 14: Plot of the effective roughness exponent C,{h,L) in 
the ordered phase, for h — 2.0 with (100) oriented faces, and 
h = 1.0,2.0,2.1,2.4, with (111) oriented faces. The samples 
have X = rL^'^ layers (of area L^ for (100) and area \/^L^ 
for (111).) For all samples in the ordered phase, the exponent 
approaches ^ = 0.66 ± 0.03 as L ^ oo, consistent with the 
expected ^ — 2/3. For comparison, data for the disordered 
phase is included; the apparent exponent decreases for large 
systems when h > he- 



2/3 

where the Xi^2 are chosen to have the values rL^ 2 ,with 
r fixed at close to unity. Assuming that C is indeed 
near 2/3, this choice ensures that a typical wall is found, 
rather than the best of a set of ^ L^~'^ possibilities that 
would result from using a sample that was much longer 
than L^ in the x direction. Such a sample shape would 
result in the same asymptotic value for C, but would have 
(probably lo gar ithmic) corrections to scaling. As can be 
seen in Fig. ^, the effective exponent appears to con- 
verge to C = 0.66 ± 0.03 in both geometries. Note that 
even with the appropriate anisotropic scaling, the correc- 
tions to scaling are large for samples up to L = 100 with 
the corresponding X ~ 20. 
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FIG. 15: Estimates di(h,L) of the box-counting fractal di- 
mension di — di{hc, 00), for the (100) orientation of con- 
trolled faces. Comparisons of the H — configuration are made 

with the and +-|- configurations in boxes of volume B^. 

If the -I — configuration differs from both of the others, that 
box is considered part of the domain wall. The finite loga- 
rithmic derivative of the scaling of the number of such boxes 
with sample size L yields the estimates shown with the lines 
connecting data points with the same h. The error bars repre- 
sent la statistical uncertainties. For h = 2.27, the dimension 
estimate converges to dj = 2.24 ±0.03, the error being a com- 
bination of statistical error and systematic errors (« 0.02) 
caused by finite size effects and uncertainties in he- 



D. Incongruence box-counting interface exponent 

di 

For an alternative measurement of the dimension of the 
domain walls at criticality, we have used a box counting 
method. In this method, we compare the configuration 
given -I — (or — +) boundary conditions with both ++ 

and configurations. This is done at various scales w, 

by partitioning the sample into (L/B)^ cubes of volume 
B^. If the configuration with twisted boundary condi- 
tions differs from both -\ — \- and in a given volume 

B^, that cube must intersect the domain wall. But this 
wall will not include any boundary of frozen regions that 
is isolated from other broken bonds by a distance of at 
least B. In particular, when S = 0, the number of such 
intersecting boxes N{B, L, h) will be zero for B smaller 
than the size of the frozen interior region. For exam- 
ple, the -I — and — h configurations in Fig. Wn are locally 



congruent everywhere with either the or the H — \- con- 
figuration. Thus only for boxes larger than the width of 
the interior region will a domain wall be apparent. 

The scaling of the number of intersecting boxes 
N{B, L, h) with L gives an alternate estimate of an ef- 
fective fractal dimension which we call dj{h,L), antici- 
pating that N ^ (L/B)'^' at the critical point (see Fig. 
pT](b).) Using the same form of the discrete logarithmic 
derivative between scales L and 2L as in Eq. (|3^) gives 
the effective exponent di{h,L), as summarized in Figs. 
M and n^. This estimate yields a constant at large L, 
within statistical errors, for h = 2.27 ~ he and gives a 
value dj = 2.24 ± 0.03. 

We note that a useful compatible definition for dj can 
be based on bonds rather than spin blocks: count the 
number of bonds that are broken with the H — or — h 
BCs that are satisfied with both the -|--|- and BCs. 
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The number of such bonds TV' should have the same scal- 
ing form as N does for fixed B. We have used this bond 
definition in a smaller number of samples and find results 
for dj(h, L) at large L consistent with the spin block def- 
inition of di defined above. 



E. Exchange stiffness exponent dj 

A third measure that we have used to study domain 
wall geometry is the contribution of the exchange energy 
to the stiffness S. This we denote Sj. It is the signed 
sum of the broken bond weights, counted as negative for 

the and -I — h configurations and positive for the — h 

and H — configurations. 

As in computing E, using the symmetrized energy dif- 
ferences reduces boundary effects. If E = 0, then Sj = 0, 
for example, though comparing the configurations with 
H — h and H — boundary conditions in such a sample will 
reveal a domain wall while comparing those with -|--|- 
and — h boundary conditions will reveal a second en- 
tirely distinct domain wall. Either of these domain walls, 
along with a portion of the frozen spins that make up the 
boundary, would be counted in the method which yielded 
ds- But in this symmetrized measure from Sj, the signed 
sums would cancel, so that neither domain wall would be 
counted. Similarly, when the box size B is smaller than 
the size of the frozen interior, the measure used to find di 
would also be zero. Note, however, that Ej does include 
some of the boundaries of the frozen regions but it does 
so with signs that can be either positive or negative. In 
the ordered phase, then, the exchange stiffness Sj will 
include contributions from the region between the two do- 
main walls that occur, contributions that would not have 
been included in the other methods. The three proposed 
measures are thus potentially all different, especially off- 
critical, but perhaps also at criticality. 

At the critical point the exchange energy part of the 
symmetrized stiffness will have contributions from the 
domain walls with holes, ~ ijd.s-(i/y ^ equivalent to the 
box counting measure of the domain wall, as well as con- 
tributions from parts of the boundaries of the frozen re- 
gions. The simplest expectation is that the the contribu- 
tions from the frozen region boundaries will be random 
in sign and thus less important in toto. 

The mean of T,j{h, L) can be used to compute a fractal- 
dimension- like quantity, dj, for the interface via the as- 
sumption that Sj ~ L'^-' at he- The scale dependent 
effective exponents from our data at h = 2.27, shown in 
Fig. 16, yield an estimate of dj = 2.18±0.02 that appears 
to be slightly smaller than the other two dimensions dg 
and di. 

One advantage of the exchange energy is that we can 
relate this measure of the fractal dimension of the domain 
walls at the critical point to the other exponents. If a 
small additional exchange SJ is added to the Hamiltonian 
(or equivalently if all the random fields were decreased in 
magnitude by a uniform small amount) then the change 
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FIG. 16: Estimates dj{h,L) of the scaling of the exchange 
contribution to the stiffness defined as the total signed surface 

area of the changes between -I--I-, , — \- and H — boundary 

conditions. The logarithmic derivative of T,j{L) gives the 
values shown with the lines connecting data points with the 
same h. The error bars represent la statistical uncertainties. 
For h = 2.27, the dimension estimate converges to dj = 2.18± 
0.03. 
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FIG. 17: Estimates d,{h = 2.27, L), di{h = 2.27, L), and 
dj{h = 2.27, L) of the fractal dimensions using controlled 
boundary surfaces in the (111) plane of the cubic lattice which 
are rhombi with sides of L spins. The number of layers in the 
sample, including the boundary planes, is L or 4L, as shown in 
the key. The error bars represent la statistical uncertainties 
with the lines connecting data points with the same h. The 
dimension estimates converge to ds = 2.30 ±0.02, dj = 2.25 ± 
0.05, dj — 2. 18 ±0.03; these are consistent within errors with 
those from Figs. |l2[ |l|, and |I|. 



in the stiffness would be simply 



6Y. 



SJ_^ 
J ' 



(35) 



Since S - L^ while Sj ~ L'^' with dj > 6, the 
change in the stiffness will become of order the stiff- 
ness itself and thus strongly modify the system when 
L ~ (SJ)~^^^'^-'~^\ This crossover length is thus a mea- 
sure of the correlation length, ^, and we thus expect the 
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exponent equality 



— = rfj — I 

V 



(36) 



This can be derived directly from the scaling form Eq. 
(Eq) by differentiating with respect to J (equivalently 
with respect to —K) and noting the thermodynamic iden- 
tity between derivatives with respect to coefficients of 
terms in the Hamiltonian and expectations of the corre- 
sponding term. (Note that this is closely analogous to 
the relation between v and the energetic part of the in- 
terfacial free energy at conventional finite temperature 
critical points.) Assuming the scaling relation Eq. (pq) 
would give v = 1.45 ± 0.10, a slightly different, but con- 
sistent, value of v than that from the scaling of the total 
symmetrized stiffness S. 
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FIG. 18: Mean absolute magnetization per spin, \m\, plotted 
vs. h for various L, for periodic boundary conditions. 



F. Comparison of domain wall exponents 

Due to the subtleties introduced by frozen islands and 
the representation of the Hamiltonian as the sum of do- 
main wall and random field components, there are three 
natural measurements of the domain wall surface and the 
domain wall contributions to the stiffness. Each measure 
has its own physical meaning. We will argue in Sec. pXB 
that the difference between d^ and dj is due to frozen 
islands, and hence ds — dj should be related to (3/i^. 



VII. MAGNETIZATION 

Having established the location and order of the tran- 
sition, we now focus on an apparently problematic quan- 
tity: the magnetization. Given sets of ground state spin 
configurations {si}, the distribution of magnetizations 
can be studied as a function of L and h. In order to bet- 
ter understand the large volume limit, we have computed 
the magnetization distributions for five different bound- 
ary conditions: all boundary spins fixed to a single value, 
either all positive or all negative (F); boundary spins 
fixed at independent random values (i?); open boundary 
conditions (O); periodic boundary conditions (P); and 
a combination (Q) with conditions P, O, and R along 
each of the three axes . The fixed spin boundary condi- 
tions will tend to favor ferromagnetism, the random will 
tend to favor a disordered phase and the combination Q 
appears to significantly reduce some finite size effects. 

We first describe our results for the mean of the abso- 
lute value of the magnetization density 



= iE^''i^' 



(37) 



for cubic samples with periodic boundary conditions (P). 
Fig. pSl is a plot of our data as a function of /i, for vari- 
ous L; the magnetization drops off quite steeply near h^. 
Fig. n9 shows the magnetization as a function of system 



size, along with its discrete logarithmic derivative, which 
yields an effective scale dependent exponent. To within 
errors, the magnetization is consistent with power law 
scaling. 



L-P'\ 



(38) 



with P/v = 0.012 ± 0.004. For h < he, the magnetiza- 
tion appears to approach a constant (e.g., to(2.255, L — > 
oo) « 0.952). For h > he, the effective exponent de- 
creases significantly as L increases. 

For further analysis, we characterize the distributions 
of m by the average over samples of the square of the 
magnetization per spin, m?, and the root- mean-square 
sample-to-sample variations of the square of the magne- 
tization 



(m 



2^2 



(39) 



Our results for A„j2 are shown in Fig. |2^. As L is in- 
creased, the peak magnitude of A^^ is seen to decrease 
for some boundary conditions F, O and P, while it in- 
creases for others R and Q. For boundary conditions R, 
O, P, and Q, the peak heights appear to be converging 
to a similar fixed value, bracketed from above and below 
by the different sets of data. In general, we would expect 
that the height of these peaks would scale for asymptot- 
ically large sizes as L~'^^/'^; the data are thus consistent 
with either /3 = or with a very small (3. 

One can estimate the location of the transition by fit- 
ting the data for A„2 near the peaks at five or more 
values of /i to a Gaussian form. (The Gaussian gives a 
better fit than a parabolic form over a larger range of 
h, though either form should give the same limit for h 
near enough to the peak and L large enough.) The fit- 
ted location of the peaks is extrapolated for all boundary 
conditions as a function of L. We obtain agreement of 
the extrapolations for 1.3 < v < 1.45 with a value of 
he — 2.272 ± 0.004, consistent with the value from Pq 
and other estimates. We believe that this independent 
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FIG. 19: (a) Mean absolute magnetization per spin, \m\, 
plotted vs. L for various h (with periodic boundary con- 
ditions.) The solid line is |m| = (1.0009)L~""". (b) 
The discrete logarithmic derivative Aln(|7n(L)j)/Aln(L) = 
\n{\m{L')\/\m{L))/ln{L' /L), vs. VlU, with L' ^ 2L. This 
is used to directly estimate P/u, yielding /3/i/ — 0.012 ±0.004, 
where the error bars are dominated by the range of values for 
he obtained by fitting over sizes up to L = 256. 



estimate is relatively precise and robust, due to the va- 
riety of boundary conditions used, with the variation in 
the results giving an estimate of systematic uncertainties. 

For the fixed spin boundary conditions, F, the peak 
magnitude of m^ is apparently converging to a different 
value (note that the magnetization near the surfaces will 
vary less than with the other boundary conditions.) If 
either these data, -F, or the periodic boundary condition 
data, P, at the critical point are used (rather than the 
data near the peak) then a smaller value of A„2 is found, 
roughly the same although apparently still distinct for 
these two cases. 

Collectively, our magnetization data would appear to 
suggest a picture of the transition that is consistent with 
that of reference Ref.[l^: three possible "states" at the 
critical point, "up", "down", and "disordered" as would 
occur at a first order para- to ferro- magnetic transition. 
As we shall see, however, our other data and further 
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FIG. 20: Magnitude of sample-to-sample fiuctuations A^2 
in the mean square magnetization per spin, as a function of 
random field strength h, for various system sizes, L, (a) for 
fixed Si = -1-1 boundary conditions (F), (b) for open (free) 
boundary conditions (O), (c) for random fixed spin bound- 
ary conditions (i?), (d) for periodic boundary conditions (P), 
and (e) for mixed periodic, random fixed, and fixed boundary 
conditions, one along each axis (Q). The curves are fit lo- 
cally with Gaussians in the regime where A„2 is greater than 
approximately 3/4 of its peak value. Extrapolating the peak 
locations to L = oo gives a best fit value oiv = 1.38±0.08 and 
he — 2.272 ± 0.004, with the dominant errors being system- 
atic errors arising from variations in the extrapolated values, 
presumably due to corrections to scaling. The lines shown are 
spline fits to visually organize the data. 



19 



thought suggest that this picture, while a very good ap- 
proximation, is not correct. We will argue that in fact 
(3 is small but non-zero and thus in astronomically large 
samples the magnetization will decay slowly to zero at 
the critical point but with the scaling functions for the 
distribution of the magnetization (and their moments) 
depending on the type of boundary conditions as is tha 
case for pure systems at conventional critical points.o 
Data suggesting this is presented in the next section. 



VIII. SPIN CLUSTERS AT CRITICALITY 

The distribution of the magnetization studied above 
gives some information about the ground state correla- 
tions of the RFIM. But because ground state correla- 
tions between Ising spins are controlled by the proba- 
bility that a pair of spins of interest are in opposite di- 
rections, the observation that the magnetization at the 
critical point tends to be rather close to saturation sug- 
gests that the loss of correlations as the random field is 
increased through the critical point may be associated 
with rather rare events. In this section, we investigate 
the nature of the effect that we believe gives the domi- 
nant contribution: the occurrence of connected clusters 
of spins of the one sign completely surrounded by spins of 
the opposite sign. Because all of the exchanges are ferro- 
magnetic, such isolated inverted clusters will, a fortiori, 
not change when the boundary conditions are inverted; 
either the spins surrounding them will flip in which case 
they will be content the way they were, or the surround- 
ing spins will not flip and the spins in the cluster will 
be isolated from the boundary condition change. Thus 
these isolated spin clusters are frozen. 

We have computed the statistics of the domain walls 
that enclose isolated spin clusters in 5000 or more sam- 
ples of system sizes up to 128'^ and 1000 samples of size 
256^ at /i = 2.27 



shown in Fig. 21 



V he- A slice of a configuration is 
Statistical errors in the dimension 
estimates and number distribution were computed by 
a bootstrap method (resapipling the statistics over the 
computed configurations) ;c2l the error bars indicate the 
estimated rms fluctuations in the statistics at each clus- 
ter size. 

We pote that previous work by Esser, Nowak and 
Usadelc3 studied the domain structure for a single sam- 
ple size. They address questions of percolation in the 
3D Gaussian RFIM, but they claim that the cluster dis- 
tribution is not broad. We find, in contrast, that there 
is a broad tail, which, though weak for smaller systems, 
becomes more important as L increases at the critical 
point. To directly contrast with the results of Ref.^, we 
find that the sum of the volume fraction of the two largest 
clusters, though near 1, slowly decreases as L increases, 
at /i = 2.27 K, he- The transition separates a state with 
one infinite connected set of spins of the same sign from 
a disordered state with two antiparallel incipient infinite 
clusters. 
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FIG. 21; A slice of a spin configuration in a 256 sample at 
h — 2.27. The dark squares indicate an up spin. The nesting 
of spin clusters can be seen here — the number of levels of 
nesting, the "depth", of the full configuration is fc = 3. The 
domain walls are determined by working recursively inwards 
from the majority (down) spin cluster. The surface of each 
cluster is taken to be the outer surface and does not include 
the surface of subclusters. 



A. Cluster surface 

For each cluster, the total volume, v, — which includes 
the volume of "holes" of opposite spin — was computed, 
as was the surface area, a, of the cluster; the number 
of unsatisfied nearest neighbor bonds that separate the 
cluster from its surrounding region of opposite spin. The 
domain walls are found recursively, taking as the initial 
surrounding region the majority spin cluster, which typ^ 
ically occupies > 97% of the volume at he for L < 256. ej 
Binning the clusters by volume v, logarithmically spaced 
by powers of 2, averaging the surface area in each bin, 
and taking the discrete logarithmic derivative gives an 
estimate of the fractal dimension of the cluster surfaces, 
dimension d'^{h, L, v). As indicated in Fig. ^, at the crit- 
ical point the surface area appears to scale as -;j0-755±0.07 
for intermediate size clusters with 1 -^ v ^ L^. The 
error in this exponent includes both statistical error and 
the apparent uncertainty of corrections to scaling that 
are affecting the convergence to a constant value. This 
value is little affected by the estimate of the location of 
he (varying h changes the number of clusters, but, within 
the uncertainty of he, does not affect the geometry of the 
domain walls.) We have verified that the volume enclosed 
by the domain walls separating opposite spins scales in 
a manner numerically consistent with this volume being 
nonfractal: 
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FIG. 22: Dependence of the surface area (number of broken 
bonds) of cluster boundaries on the enclosed volume, v, ex- 
pressed as an effective exponent ds(/i, L, v), at /i = 2.270 ~ he, 
for L — 32, 64, 128, 256. The cluster surface area scales as 
^o.755±o.oo7 £qj. ^j^g largest clusters that are not affected by 
finite size effects, yielding a fractal dimension d^ = 2.27(2) 
for the cluster surfaces. 



with £ being either the geometric mean or the maximum 
of the lengths of the sides of the minimal rectilinear box 
that encloses the cluster. The extrapolation of d'^{h, L, v) 
to large L and £ is therefore consistent with clusters hav- 
ing typical diameter £ ~ v^''^ and typical surface area 



with 



2.27 ±0.02 



(41) 



(42) 



a fractal surface dimension consistent within the statisti- 
cal uncertainties with our estimates of the fractal dimen- 
sions ds and dj of the domain walls induced by changing 
boundary conditions at the critical point. In particular, 
this surface dimension bears a close resemblance to the 
dimension of the spanning surface which we denoted ds ; 
thus we conjecture that 



dl = ds 



(43) 



FIG. 23: Fraction of the volume p{v) occupied by clusters 
of volume between v and ev, ai h = 2.270 ~ he, found by 
normalizing the data binned according to powers of 2 (i.e., 
dividing the volume fractions in the [v, 2v] bins by In 2.) The 
solid line is p{v) = 0.0019 -I- (0.0017)^;"°-^^, one of the trial 
fits used to extrapolate to large v. 
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FIG. 24: Fraction of the volume p{v) occupied by clusters 
that are contained in rectilinear volumes ("boxes") between 
V and ev, at ft = 2.270 ~ he, found by normalizing the data 
binned according to powers of 2 (i.e., dividing the volume 
fractions in the [w,2w] bins by In 2.) For large v, p{v) —> 
0.0019 ± 0.0002, if he = 2.270. 



B. Cluster density 

New information is given by the densities of the clus- 
ters as a function of their size, in particular their dimcn- 
sionless volume fraction 



p{v) 



-Prob[site e cluster of size in {v, v + Sv)]. (44) 



From the data in Fig. E3^ (and for the slightly different 
measure of Fig. 24, where v is the volume of the smallest 
parallelpiped of fixed orientation enclosing the cluster) 
we see that clusters that are neither too small nor lim- 
ited by finite size effects — roughly a decade in length 



scale for L = 256 — occupy an approximately scale in- 
dependent volume fraction. A comparison of the cluster 
distributions for nominally off-critical values of h, as seen 
in Fig. p3, shows how p{v) depends on h. From these 
plots we infer a large-u limit of 



p{v) 



Po 



0.0019 ±0.0004. 



(45) 



We cannot, of course, rule out a slow decrease of p{v) to 
zero for large volumes, especially as our effective range of 
length scales here is less than for other quantities because 
of the restrictions due to finite size effects. But we can 
understand on the basis of our other observations why 
one should expect a small but non-zero value for pa^- 
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FIG. 25: Fraction of the volume p{v) occupied by clusters of 
volume between v and ev, found by normalizing as in Fig. 
G3. Here, the volume fractions are plotted for L — 32, 64, 128, 
with h = 2.255, and L = 256, with h = 2.255, 2.270 and 
2.285, to indicate some of the effect of changing L or ft on p{v). 
The data for h = 2.255 apparently converges at large L to a 
well defined distribution that has a finite-u cutoff, consistent 
with a finite correlation length in the ordered phase. This is 
in contrast with the data for h = 2.27, the putative critical 
point, for L < 256, as seen in Fig. ESl For h = 2.285, in 
the disordered phase, large volume clusters start to occupy a 
larger fraction of the volume than smaller clusters, for L — 
256. 



IX. SCALING 

In this section we pull together our various results 
about domain walls, stiffness, magnetization, and in- 
verted spin clusters and show how they are all consis- 
tent with a simple picture of scaling behavior at a zero 
temperature phase transition. 



A. Critical correlations 

At the critical point, the energy cost of domain walls is 
typically sufficiently large that almost all cubical samples 
— about 96% of them — would rather have no spanning 
(or other large scale) domain walls unless forced to by 
boundary conditions. But in a small fraction of the cu- 
bical samples the random fields in the central region are 
sufficiently strong that they force the system to have two 
domain walls for one of the two "ferromagnetic" {++ or 

) choices of boundary conditions. In samples that 

are twice as long, this occurs much more frequently as 
evidenced by the increase, on going from cubical to elon- 
gated samples, in the probability, Pq, that the stiffness 
vanishes. Although whether such a pair of walls is fa- 
vorable generally depends on both the random fields in 
the whole system and the local behavior near the walls, a 
crude picture of what is going on can be drawn by assum- 
ing that the wall energies are relatively local and weakly 
dependent on each other. We restrict consideration for 



now to the critical point. 

First consider a system of dimensions ^Lx Lx L with 
the boundary conditions imposed on the faces perpendic- 
ular to the short axis. Assume that the probability that 
the single wall energy, E^y^ = E^ —Ej^+, in such a sys- 
tem is negative is g <C 1. Crudely, for two walls to be fa- 
vorable in a cubic system with -| — \- boundary conditions, 

<0 
<0 



as is needed to make E = , one must have both E^ 
for the left half of the system and E^ = E ^ — E- 



w 



for the right half of the system. Naively, this occurs with 
probability of order q^. (More precisely, one of the Ew's 
could be positive but not by enough to dominate the 
other one; this will not change things much as long as 
the bulk of the distribution of the Ew's is skewed sub- 
stantially to the positive side of zero.) But in a system of 
length 2L rather than L, there are many more possibili- 
ties: if we divide the system into four sections of length 
L/2, one could have, for example, the second from the left 
having E^ < and the rightmost having E^^ < with 
the wall energies of the other sections being positive. As 
there are six such choices among the four sections of the 
elongated system, we expect that the chances of having 
E = will be about six times as large as in the cubical 
system — obviously a very crude approximation, but one 
that yields roughly the measured magnitude of the ratio 
Po{2L xLx L)/Po{L xLxL). Note that this picture im- 
plies that for systems that are much longer than they are 
wide, the typical number of domain walls in the ground 
state will grow linearly with the length. The roughly 
random spacing between them will lead to exponential 
decay of the end-to-end correlations in such a system, 
with a characteristic length proportional to the linear di- 
mension of the cross-section as should be expected on 
general finite size scaling grounds. 

At conventional critical points in two dimensions, con- 
formal invariance relates the exponential decay of correla- 
tions in long tubes to the power law decay of correlations 
in the bulk in infinite systems: the exponent rj is simply 
proportional to the ratio pf the width-dependent correla- 
tion length to the width.Cj In our case, there is no such 
exact relation, but one can make a qualitative argument 
that suggests a similar result. Consider a region of diam- 
eter of order £ centered on some chosen spin in the bulk 
of the sample and assume that outside of this region, the 
spins in the vicinity are -1-1. The only way that the spins 
inside the region of interest can be —1 is if there is a do- 
main wall relative to the pure "up" configuration which 
surrounds this region and has negative energy. Roughly 
speaking, such a closed domain wall must be made up 
of four or more sections which are joined together with 
each having negative (or close to zero) energy. Since the 
amount of freedom perpendicular to the area of each of 
these will be somewhat less than their linear dimensions, 
a crude approximation is that the probability of finding 
each such section is of order q and the probability of find- 
ing the total domain wall energy negative is of order g^. 
With Po[cube] ~ 0.04 ~ g^, this suggests that the prob- 
ability for finding such a spin flipped region will be of 
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order po 



0.002, in the same range as that found. 



Although the above argument is very crude — factors of 
two's or pi's could easily have been fudged in! — it nev- 
ertheless provides a suggestive connection between the 
smallness of various quantities. Indeed, it actually pro- 
vides more: a precise method for estimating [3/ v. 

Consider a single spin in the center of a large system 
at the critical point with, for simplicity, -I- boundary con- 
ditions . For each factor of e^/'^ in length scale, £, there 
is a probability Poo that the spin is an element of a clus- 
ter, flipped with respect to its surroundings, with vol- 
ume in the associated range P^ < v < el'^. There will, 
of course, be correlations between the probabilities of oc- 
currences of such inverted regions that are similar in size 
and nearby to one another. But, because inverted re- 
gions are so rare, the effects of these correlations will be 
negligible and we can assume that each range of i around 
the chosen spin is independent. A simple picture of the 
behavior then emerges: the spin of interest will be in 
a cluster of one orientation of diameter ii, which itself 
will be in a much larger cluster of the opposite orien- 
tation of size f.2, etc. with the successive sizes growing 
approximately geometrically — more precisely as a Pois- 
son process in ln(i?) with density dpoo (with d = 3). The 
probability, p||, that the spin has the same orientation as 
the largest cluster — the system size L — can readily 
be computed from the properties of the Poisson process; 
this result yields the mean value of the spin given fixed 
(+) boundary conditions 



o( + ) 



2p|| - 1 



1 



LP/^ 



with the exponent 



2dpoo ~ 0.011 ±0.003 



(46) 



(47) 



A similar argument for two spins a distance ja; — y| apart 
in an infinite system gives 



I 



|2,_y|d-2 + i) 



(48) 



with the modified "anomalous dimension" exponent for 
these untruncated zero-temperature correlations given by 



d-2 + fj = 213 /v w 4rfpo 



(49) 



This picture of droplets within droplets strongly sug- 
gests that neither spin species will percolate at the critical 
point. This is, a priori, rather surprising as the percola- 
tion concentration for a three dimensional cubic lattice 
is substantially less than half and so one might have ex- 
pected both spin species to percolate even somewhat into 
the ordered phase. The fact that they do not in this sys- 
tem is associated with the smallness of /3 and the nature 
of the critical point. 

In practice, unfortunately, the value of P/v is so small 
for the 3D RFIM that the effects discussed above will 



be ail-but unobservable even if experiments could reach 
equilibrium. But in higher dimensions, four or five, they 
might be observable numerically as relatively large sys- 
tems sizes (e.g., more than 32^) can be explored. 



B. Fractal dimensions of domain walls 

The picture developed above suggests that the vari- 
ous fractal dimensions of interfaces or domain walls at 
the critical point will not be the same but might, nev- 
ertheless, be related to the other exponents. The fractal 
dimension of the spanning interface dg (and the dimen- 
sion of the surfaces of clusters) is the dimension of a true 
surface, one with no holes in it. Such a surface cuts across 
the whole system but the sets of sites it is separating can- 
not really be thought of as belonging to different states 
— the "up" and "down" states — since, in an asymp- 
totically large system, most of the sites will be frozen 
and unaffected by the boundary conditions. In contrast, 
the incongruence box counting dimension, d/, is sensi- 
tive only to those parts of the system that are affected 
by boundary conditions: a fraction of order XjL^^^ . A 
natural conjecture is that the box counting dimension is 
the same as that of the intersection of a typical fractal 
spanning surface with a fractal set of dimension d — P/v 
yielding: 



P/u. 



(50) 



This picture is somewhat analogous to what would occur 
right at percolation in a diluted ferromagnetic Ising sys- 
tem at zero temperature: in a finite fraction of the sam- 
ples, forcing a domain wall by changing the spin bound- 
ary conditions would cost no energy, while in the rest it 
would cost an energy proportional to the area of an in- 
terface that only cuts across the fractal incipient infinite 
cluster; this interface would have dimension analogous to 
our dj. 

The exchange energy dimension, dj is sensitive to both 
the frozen and the unfrozen regions. But a reasonable 
guess is that this is dominated by the unfrozen regions as 
the frozen regions contribute random signs. This suggests 
that 



dj = di = 9 + Ijv 



(51) 



which, if correct, implies that a relation obtains between 
ds and the other exponents : 



ds = 



1 + /3 



(52) 



We should note, however, that these conjectures are 
difficult to test in three dimensions, due to the smallness 
of /3. Our estimated exponents are in slight disagreement 
with these conjectures, but corrections to scaling that are 
not apparent can be important at this level of accuracy. 
Nevertheless, our apparent values for d^ and d^ do appear 
to be larger than di and dj. 
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In higher dimensions, testing the conjectured scaHng 
relations between these dimensions and the other expo- 
nents might be feasible. It is of course possible, however, 
that further analytic understanding would imply that at 
least some of the fractal dimensions could be independent 
exponents. 



C. Specific heat 

The speeific heat of the RFIM can be experimentally 
measuredED and is of theoretical importance. Monte 
Carlo methods atJpijbe temperature have been used to 
estimate its value rwJ\ In addition, Hartmann and Young 
(HY) have recentlyEj computed the exponent a describ- 
ing the divergence of the specific heat, using ground state 
configurations. They find a value of a = —0.63 ± 0.07. 
Using the same thermodynamic assumptions, but differ- 
ent analysis methods, we find a — —0.01 ± 0.09. 

One expectsEa that the finite temperature definition 
of the specific heat can be extended to zero tempera- 
ture, with the second derivative of {E) with respect to 
temperature being replaced by the second derivative of 
the ground state energy density Egs with respect to h, 
or equivalently up to constants, J. The first derivative 
dEgs/dJ is just the average number of unsatisfied bonds 
per spin, Ej = i~'^X]<ii> ^i^r Hartmann and Young 
(HY) calculate the needed second derivative by finite dif- 
ferences of Ej{h) for values of h near he- [Ej is not ex- 
plicitly dependent on J, but changes discontinuously in 
a finite sample when the spin configuration {s.;} changes; 
the second derivative is thus a set of (5-functions which 
are smoothed by the finite differencing.) The finite-size 
scaling form assumed is that the singular part of the spe- 
cific heat C, behaves as 
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FIG. 26; A plot of Ej, the bond part of the energy density, 
for h = 2.27, as a function of L. The fit shown is of the form 
-Ej = ci - caL^""^^/"", with ci = 0.14632, ca = 0.29098, and 
{a-l)/v = -0.82. The residuals (inset) give x^ = 0.65. The 
statistical error in (a — l)/u for fixed h near he is 0.02, but 
the uncertainty in this ratio is 0.10 due to the uncertainty in 
he. 
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(53) 



HY determine a by fitting to the maximum of the peaks 
in Cs, which occur at hpea.k{L) — he ^ L^/". 

Here, we estimate a using the results for the stiffness 
from Sec. |^ and also by studying the behavior of Ej at 
he- The first estimate found by applying Eq. (pf) with 
our values of 6 and i/, is a(i) = — 0.07± 0.17. The com- 
putation from the behavior of Ej is based on integrating 
Eq. (m) up to he, which gives the dependence 



Sj,,(L,/i = M=ci+C2£("-i)/", 



(54) 



with ci and C2 constants. We have computed Ej for a 
large number of samples of various sizes and estimated 
the singular part of the sample average. We directly fit 
our data for Ej[L), at fixed h, to the form Eq. (|54|). The 
fit for the nominal he, h = 2.27, is shown in Fig. 09. 
The fitted values are (a - l)/v = -0.82 ± 0.02, where 
the quoted error is purely statistical. The fit is good for 
16 < L < 256, with x^ = 0-65 for a three parameter fit 
to five data points. This fit is also consistent with that 



FIG. 27: Plots of the discrete derivative with respect to 
ln(L) of Ej, for h = 2.255, 2.270, and 2.285. The solid lines 
show power law fits for L > 30, with slopes 1.21, 0.84, and 
0.60, respectively. Using the error estimate for he, this gives 
{a — l)/v = —0.84 ± 0.10, consistent with the results from 
plots as in Fig. Em 



found from taking the derivative of Ej with respect to 
ln(L), 



dE, 



d(ln L) 



lia-i)/^ 



(55) 



at ft, = 2.27, which removes the need to fit to ci, but 
introduces larger uncertainties, due to the derivatives. 
This data for h near he is displayed in Fig. E^. By varying 
h {h ~ 2.255, 2.280), we estimate the systematic errors, 
given our uncertainty in he, arriving at the value [a — 
l)/i/ = -0.82 ± 0.10, which, using ly = 1.37 ± 0.09, gives 
a second estimate a(2) = —0.12 ± 0.16. 

Besides the uncertainties in he, this result for a is af- 
fected by finite size corrections. We now argue that these 
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corrections can be reduced by extrapolation and that a 
connection exists between a, d^ and filv: Ej, being the 
bond part of the energy density, is simply given by the 
density of domain walls , whose scaling can be found from 
the results in Sec. VIII, Namely, taking the surface area 
of clusters to scale with linear size £ as A ^ d'^^ ^ yd.jd 
and using the constant limit for the distribution of vol- 
umes p{v)d[l'a{v)] at large v, the domain wall density in 
a finite sample is found by integrating the wall density, 
taking into account intersections between the scales, over 
L up to the system size, giving 



E 



J,s 



Tdl-d-P/i' 



(56) 



This exponent can be justified by considering the change 
in Ej,s upon doubling the system size. With finite prob- 
ability {poo In 2), an extra domain wall of scale L will be 
introduced. The connected surface of the domain wall 
will have area L », but the increase in domain wall area 
will be smaller, as the domain wall will intersect frozen re- 
gions. The fraction of the sample that is not frozen scales 
as L^^/^ at criticality; the intersection of the new wall 
and the unfrozen region therefore scales as ~ L'^^~^/'^ , so 
that the expected fraction of newly broken bonds (com- 
pared with the smaller sample) is ^ Ldi-d-p/v ^ (The 
domain wall intersects frozen regions that did not have 
surfaces, as they were embedded in like spins, adding to- 
tal area, and also intersects frozen clusters that had sur- 
face area, removing total area, but these contributions 
average to zero.) This argument implies the scaling rela- 
tion 



1 



d - 13 /v. 



(57) 



Note that this relation is consistent with modified hyper- 
scaling, Eq. (@), and the conjectured relationships among 
the domain wall dimensions, Eqs. (^) and (p^). Apply- 
ing this result to our data, we find 



{a-l)/v= -0.74 ±0.02, 
giving our best estimate 

a = -0.01 ±0.09. 



(58) 



(59) 



Note that the magnitude of the (3/v contribution is small 
compared with the error. 

Our result for a is in marked disagreement with the 
value from HY. The scaling assumptions for our and HY's 
analysis are identical. It may be that one set of results is 
more strongly affected by finite size errors, though we do 
fit larger values of L. We note that the value of a that 
we find using Ej is extremely sensitive to the assumed 
value of he and that the uncertainty in h^ dominates the 
error estimate. A change of he by She — 0.01 gives a 
change 5[{a — l)/v] ~ 0.2 or 5a « 0.3. We are fitting 
for values near h^^ whereas the peaks in C found by nu- 
merical differentiation are somewhat above he- In Sec. 
M, it is found that the convergence to a scaling function 



for h — he more than a couple of times L~^l^ ^ where the 
peaks in C are, is slow compared with the convergence 
at ft, = he- 

We use here two independent data sets to arrive at 
our estimates for a: (a) total stiffness measurements on 
isotropic and anisotropic samples, with fixed BC's on two 
walls, applying finite size scaling, and (b) the measure- 
ments of the bond part of the total energy, _Ej, using 
periodic isotropic samples, and fitting using a finite size 
scaling form. For (b), we analyze the samples in two 
ways: directly extracting Ej data and also estimating the 
asymptotic scaling using the d% and (3 jv measurements. 
This latter method is least sensitive to uncertainties in 
he. 



D. Deviations from criticality 

As the system is taken away from the critical point, 
the nature of the spin clusters and correlations changes 
in a straightforward way. If the exchange is increased, 
driving the system into the ordered phase, then the large 
inverted droplets, which typically have gained energy of 
order i at the critical point, will usually have this en- 
ergy gain overcome by the extra exchange energy cost of 
order liJl'^j when i is greater than the correlation length, 
^. Large inverted regions will be exponentially rare on 
length scales longer than ^. 

If the random field is increased or the exchange de- 
creased to drive the system into the disordered phase, 
we can no longer simply focus on the inverted regions 
that already exist at the critical point but must also con- 
sider putative inverted regions that could exist. In any 
region with diameter of order i, there will be, at the crit- 
ical point, an excitation that flips of order ^'^ spins for 
a typical energy cost of order £ (more precisely it will 
only flip of order £'^~P/'^ because of the frozen regions 
within it which are not sensitive to the boundary of the 
region). Since decreasing J will decrease the energy cost 
of this excitation by an amount of order SJL'''' , a good 
fraction of these "excitations" will have negative energy 
and thus occur spontaneously at scales of order ^. On 
this and larger scales, the orientation of the spins will be 
determined primarily by the local random flelds within a 
distance of order ^ of the spins of interest. 



E. Thermal fluctuations and excitations 

The effects oLlhermal fluctuations have been dis- 
cussed elsewhereBH in the general framework of a zero- 
temperature random field critical fixed point. We will 
thus restrict ourselves here to a few comments in light of 
the present more detailed picture. 

At the critical point, as has been outlined above, there 
should be potential excitations with energy of order £^ 
around each point, an independent one for roughly each 
factor of two in length scale, £. Since the energies of 
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these are random there is a finite probability density that 
the energy of any given one of them is near zero — in- 
deed there would have been ones with negative energy 
but these give rise instead to the inverted clusters in the 
ground state. The thermal fluctuations are dominated 
by the rare active excitations whose energy is within of 
order T of zero. Because 9 is positive, the active excita- 
tions with diameters of order £ 3> 1 occupy only a small 
fraction - of order T/l^ - of the volume. But this small 
active fraction dominates the correlations, in particular 
causing the thermal fluctuations of the spin-spin correla- 
tions, the truncated correlations, to decay as 



{{Sx- {Sx)){Sy- (Sy))) 



T 



\^_y\d-2+rj 



(60) 



where the exponent rj is related to that of the zero- 
temperature correlations Eq. ( [48| ) by 



T] = fj + I 



(61) 



the extra factor oiT/\x — yf coming from the prob|ibility 
that both spins are in the same active excitation.H 

In general, except for fluctuation quantities such as 
the truncated correlations, the statements that we have 
made about zero temperature will hold with minor (if 
sometimes subtle) modifications provided one considers 
always free energies instead of energies. 

One effect which must be mentioned, however, is the 
"hypersensitivity" to changes along the critical line — 
sometimes, rather misleadingly, referred to as "chaos" . 
As long as 9 < d/2, which we believe it probably is, 
although only barely so, which spins have which orien- 
tation at the critical point will depend, on sufficiently 
large scales, |-ejeU-emely sensitively on where one is on the 
critical line.c2lcj Unfortunately, due to the smallness of 
d/2 — 9, this effect is unlikely to be observable in three 
dimensions but may be in higher dimensions for which 9 
is expected to deviate more significantly from d/2. (In 
six dimensions and above, 9 = 2.) 



X. GROUND STATES AND SENSITIVITY TO 
BOUNDARY CONDITIONS 

The simple picture of the random field Ising system ex- 
hibits two phases with a single transition between them: 
an ordered phase in which a typical spin is aligned with 
others far away; and a disordered phase in which the mag- 
netization is zero and the orientation of each spin is deter- 
mined locally by the random fields in its vicinity. In the 
ordered phase, h < he, spins have long range correlations 
and there are both "up" and "down" states, although do- 
main walls can be introduced that divide the system into 
up and down regions. In contrast, when h > he, the spin 
correlation function is short ranged, with characteristic 
scale £, ^ {h — hc)^'^ and there is only one state; because 
of the locality, large scale domain walls do not exist in 
this phase. 



But it is interesting to ask, by analogy with spin glasses 
and other systems with quenched randomness, whether 
the random field Ising system could be more complicated, 
especially near to the critical point. In order to address 
this, we must characterize the macroscopically distinct 
states in an infinite system: is there, as the simple picture 
would suggest, simply one state in the disordered phase 
and two in the ordered phase? Or is the behavior more 
subtle? 

It has been claimed in the literature that "replica sym- 
metry breaking" calculations show the existence of an in- 
termediate glassy phase, where many solutions with dis- 
tinct local magnetizations coexist for a finite range of 
parameter, between the paramagnetic and ferromagnetic 
phases.Ql3 But what does this mean? Indeed, what does 
one mean by "ground states" in an infinite system with 
random couplings? Furthermore, if one answers these 
questions, what is the connection between multiplicity of 
infinite system ground states and the notion of "replica 
symmetry breaking"? 

To consider these questions, it is simplest to restrict 
consideration to systems, such as the RFIM with Gaus- 
sian random fields, in which the finite system ground 
states for given boundary conditions are non-degenerate 
with probability one. (Otherwise one gets into the com- 
plications of ground state entropy as in dilpted antiferro- 
magnets in a field and the bimodal RFIMJ^H but these 
issues are distinct from the basic questions of "states" on 
which we focus.) 



A. Infinite system ground states 

A ground state of an infinite system with finite range 
interactions is a configuration whose energy cannot be de- 
creased by changing any finite collection of spins. Equiv- 
alently, a ground state can be thought of as the limit of 
a sequence of finite system ground states of larger and 
larger subsystems, generally with appropriately chosen 
boundary conditions on each size. Thus the set of all 
ground states for a specific infinite system, is the set p£ 
all distinct limits of sequences of boundary conditions.Ell 
For two ground states to be distinct, they must be distin- 
guishable within some finite distance of the origin: if the 
finite system ground states differ only in regions whose 
distance from the origin grows without bound as the sys- 
tem size increases, then the infinite system ground states 
are the same.Ej All infinite system ground states have the 
same energy density but comparing the energy of a pair 
of ground states is not generally well-defined. 

Many of the subtleties involved in considering infinite 
system ground states come to the fore in the ordered 
phase of the random field Ising model. If we take the limit 
of larger and larger systems with open (i.e. free) bound- 
ary conditions centered, for example, on the origin, then 
the finite system ground states will not approach a limit! 
This can be readily understood in terms of the "up" and 
"down" states which we know exist in the infinite system 
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- albeit with some finite density of misaligned spins.EJ A 
given finite sample of volume V will typically have ran- 
dom fields whose net effects are to cause an energy dif- 
ference between the up and the down states which is of 
order hvV. Thus the ground states with open bound- 
ary conditions will alternate randomly from mostly up 
to mostly down as a function of (the logarithm of the) 
system size. Of course, the up and down states can be 
found by either taking the appropriate subsequences with 
open boundary conditions, or by taking -I- or — boundary 
conditions on all sizes. The problem of energy compari- 
son is now clear: which of these two states has the lower 
energy in a specific infinite system? This is manifestly 
ill-defined, indeed, because of the effects of the bound- 
ary conditions, it is not possible to uniquely define the 
energy of an infinite system ground state to higher accu- 
racy than of order the surface area of the region under 
consideration. 

We can, however, compare the energies of some pairs 
of infinite system ground states even in random systems. 
In high dimensions, greater than three, one can make 
ground states in the ordered phase of the RFIM with a 
domain wall that passes near the origin with a chosen ori- 
entation by putting -|- boundary conditions on half of the 
boundary and — on the other half. If the random fields 
are weak enough (in four and five dimensions, or with 
arbitrary randomness in the ordered phase in six or more 
dimensions), the domain wall will be flat on large scales 
with only finite typical deviations from planarity and its 
position and orientation can then be fixed by its inter- 
section with the boundary which is forced by a "seam" 
between -I- and — areas of the boundary conditions. The 
infinite system domain wall state so constructed will be 
stable to changing any finite collection of spins, but, in 
a well-defined sense, it has higher energy than either the 
up or the down states. As the domain wall costs energy 
per unit area, if one looks at a sufficiently large region 
that overlaps the domain wall — say cubical with v = t^ 

— the difference in energy between the domain wall state 
and the up state will be of order J£'^~^ ± hi'^^^ which is 
positive almost surely in the limit of large £. 

In contrast to the higher dimensional case, in the three 
dimensional RFIM of primary interest, one cannot make 
infinite system domain wall states straightforwardly even 
in the ordered phase. If one tries to set up a domain wall 
that is, say, horizontal in a system of size L x L x L, one 
will find that the wall wanders in the vertical direction 
away from the plane determined by the boundary joint by 
a random, sample and subsystem-size dependent amount 
of order L'^ with ( — 2/3.lZIl§E^ No matter how one 
adjusts the boundary seam, one is unlikely, in the large 
system limit, to be able to force the wall to be both near 
the origin and nearly horizontal. Thus the sequence of 
domain wall forcing boundary conditions will, in the or- 
dered phase, contain one subsequence which converges to 
the up state, another which converges to the down state, 
and, almost surely, no other convergent subsequences. 
(There are subtleties, which we will not go into here, if 



one allows a wall in the ordered phase to have any con- 
figuration dependent orientation; these will be addressed 
in Ref. |l[) 

The crucial question that we would like to address here 
is whether there exists more than one infinite system 
ground state either at the critical point or slightly into 
the disordered phase. In principle, to investigate this one 
would need to study all possible sequences of boundary 
conditions, obviously an impractical task. In practice, 
one must restrict consideration to some small subset of 
boundary conditions and try to extract useful informa- 
tion about the infinite system limit by carefully studying 
the size dependence of various boundary conditions on 
regions near the origin. 



B. Numerical studies 

We have studied how the ground state configurations 
change in response both to varying the boundary condi- 
tion at fixed size and to changing the system size. We 
compare configurations for which the boundary spins are 
"open" (O), fixed positive (-I-), fixed negative (— ), and 
random fixed spins (R). For fixed size calculations, for 
each realization of the random fields we compare all possi- 
ble pairs of boundary conditions in the set {O, R, +, —}. 
We also compare ground state configurations for open 
boundary conditions on a sample of size 2L — 1 (denoted 
D) that contains a subsample of size L, with the states 
for boundary conditions O, -I- or — imposed on the sub- 
sample. (The values of L were taken to be odd for these 
comparisons, so that the origin coincides with a spin.) 
The results of all of the comparisons are characterized 
by counting how many spins differ for the two boundary 
conditions in a volume w^ centered at the origin. 

The primary emphasis of these calculations is to deter- 
mine whether changes in boundary conditions can create 
configurations that differ from those with uniform -I- or 
uniform — boundary conditions, i.e., those that produce 
the up and the down states in the ordered phase. If 
the -|- and the — boundary conditions produce identi- 
cal configurations in the deep interior, this suggests that 
there is only one state. If the probability that some other 
boundary condition produces a configuration in the in- 
terior that differs from those of both the the + and — 
boundary conditions, vanishes as L — > oo, this suggests 
that there are at most two states. 

We report here a selection of results for (a) the proba- 
bilities Po±{h,w,L) {Pji±(h,w,L)) that the boundary 
condition O (respectively, R) gives a central volume 
w'^ that differs from that for both + and — boundary 
conditions at fixed sample size L; (b) the probability 
PDo{h,w,L) that the number of differing spins within 
the window is non-zero when one compares open bound- 
ary conditions for samples of size 2L — 1 and a subsample 
of size L; and (c) the probability PD±{h, w, L) that open 
boundary conditions on the larger sample gives a cen- 
tral volume w'^ that differs from that for both -\- and — 
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boundary conditions on the smaller sample. 

The calculation of the probabilities Po± and Pji± 
(comparisons (a)) allows us to study ground states near 
he. The events of interest are those where a given bound- 
ary condition B, either B — R or B — O, gives a con- 
figuration distinct from both the -|- and the — boundary 
conditions. For h > he a.s L ^ oo, Pb± is expected to 
go to zero, since the effects of the boundary penetrate 
only a distance O(^) into the sample. For h < he, in 
contrast, most of the interior configuration is either -I- or 
— and the chances that random or open boundaries yield 
some other possibility should again decay exponentially. 
At h = he, the correlation length diverges and at this 
critical point, we expect that the probability of a domain 
wall passing near the center decays only as a power of 
the system size. Using simple . Fij:guments based on the 
fractal nature of domain walls,E3Ej the probability that 
a window of size w will intersect an object of fractal di- 
mension df scales as {w/ LY^'^^ . This is analogous to 
the probability of a domain wall in the ordered phase 
passing near the origin as discussed above. The appro- 
priate fractal dimension to use here at the critical point 
is the dimension from box-counting, dj, that we studied 
above. Basically, there is a substantial probability that 
open or random boundary conditions will, at the criti- 
cal point, induce a system spanning domain wall relative 
to the -|- and — boundary conditions. Near the critical 
point, scaling suggests the form 

PB±{h,w,L) ^ L''f-'VB[w,{h- he)L-^'''] (62) 

ior B = R or O. We plot our data for Po±, w = 3, in 
Fig. p8|and Fig. [29|(a), assuming this scaling form, taking 
he = 2.270, df = 2.25 and the best fit value ly = 1.37. 
The results for Pr±, while not shown here, are nearly 
identical, apparently converging for large L near the crit- 
ical point to an extremely similar, if not the same, scaling 
function, though the smaller L curves have slightly dif- 
ferent finite size corrections. We expect that df is equal 
to the incongruent domain wall dimension dj, as this is 
the domain wall dimension that describes changes in the 
bonds, and this expectation is consistent with our results. 

We note that taking the value df — 2.20 appears to 
give a better fit for the peak heights away from he, but 
as convergence in several quantities is poorer away from 
he, the value df — 2.25 is acceptable. Directly fitting the 
peaks for df gives a value df — 2.22 ± 0.03. 

The data for comparisons (a) should also scale with 
w for large w: V[w, {h - hejL-^^"] ^ w'^-'^fV'[{h - 
he)L^^''^]. However, we do not have enough range in 
w for w » 1 to confirm this; for w small, discreteness 
effects will prevent a collapse. For L — 129, the data do 
collapse well for w = 65, 33, 15, assuming the above scal- 
ing form and the best fit values of he, dj, and v. Note 
that similar finite w effects were also seen in the data of 
Ref. ^, where large w was needed to see convergence to 
power law behavior in w, though scaling worked well for 
fixed w with L 3> 1. 

Comparison (b) compares open boundary conditions 
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FIG. 28; Plot of the unsealed probability Po± that the 
central window of size w = 3 of a ground state configuration 
with open BCs on a given sample of size L, differs from the 
configuration in the window with both uniform + and — fixed 
boundary conditions. The lines are intended to organize the 
data visually. 
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FIG. 29; Scaling plot for the probability that the central 
window of size w = 3 of a ground state configuration differs 
from that of uniform + or — fixed boundary conditions for 
(a) open boundary conditions on the same sample of size L 
and (b) open boundary conditions on a sample of size 2L — 1. 
The values used for scaling are he — 2.270, u — 1.37, and 
df — 2.25. The probabilities scale very well near h = he, but 
the peak heights, at /i > he, converge slowly. 



on two samples of different sizes, the smaller being a 
subsample of the larger centered at the origin. In the 
disordered phase, with the local spin configuration de- 
termined by the random fields nearby, doubling the size 
of the system is not expected to change the configu- 
ration in small windows near the origin, for w <^ L 
and L ':S> i, ^ {h — he)~'^ ■ But for h < he and 
L ^ ^ ^ {he — h)^'^ , the spin orientation is determined 
by the sign of the total (effective) random field which will 
depend stochastically on the system size as discussed in 
the previous subsection. For a fully magnetized system, 
{\m\ = 1), these simple expectations yield Pdo -^ for 
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L — > oo with w fixed for h > hr, and 



Pi 



DO 



(Ttt^) 



-1/2 



da; 



dye 



(-aV2-(y+2:)Vi4) 



^ (63) 

for /i < /ic as L — !■ oo. The integral in Eq. (63) is the 
probabihty that the total random field in the volume 
(2L)^ exceeds in magnitude and is opposite in sign to the 
total random field in a subvolume L^, assuming Gaussian 
distribution of the field on length scale L with variance 
L^ . This integral gives a value Pdo = 0.384973 . . . for 
h < he. The results of our ground state studies, dis- 
played in Fig. ^, appear to be consistent with this limit, 
for h < he- This confirms the expectation that as the infi- 
nite volume limit is taken in the ordered phase, the spins 
in a fixed volume flip between two distinct conflgurations 
infinitely often — typically every factor of three or so in 
length scale. Near the critical point, the probability of 
differences in the window between the full and the sub- 
sample will be modified since with |m| < 1, there is a 
non-zero probability that the window will be contained 
in a frozen spin clusters that is unaffected by the overall 
majority random field. But this ui-dependent difference 
only becomes important near he, as /3 is so small. 

Right at the critical point the effects of frozen clusters 
on all scales should in principle suppress Pdo to zero in 
the large L limit for all w; but as it will decay only as 
l/L^/"^, this effect is hard to see. In the disordered phase, 
our data is consistent with Pdo vanishing exponentially 
for L > e 

Comparison (c) allows us to address nearly the same 
question as (a), but more directly checks that increas- 
ing the volume of the system has the effect of setting an 
effective boundary condition of -I- or — on the central 
region. The scaling collapse, shown in Fig. ^(b) is ac- 
ceptable, as it is in (a), with a scaling function similar 
to, but distinct from that for comparison (a). The re- 
sults for Pd± show that, except for a region near he that 
shrinks and decreases in probability with increasing L, 
the configuration given by the larger system with open 
boundary conditions does not produce a distinct interior 
volume from that found by imposing -f or — boundary 
conditions on the smaller system of size L. 

Taken together, these results are consistent with the 
expectations from the simple scenario for the structure 
of the states given above; there do not seem to be any 
indications of stranger behavior. Thus, in the absence 
of concrete testable predictions from those who believe 
there should be more than just the simple set of states, 
we can do no more than conclude that if they can indeed 
occur, it must be only under very subtle conditions. 
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FIG. 30: Probability Pdo that the central window of size 
(a) u) = 1 and (b) lo = 3 of a ground state configuration in a 
subsample of size L differs from that in a sample of size 21/ — 1, 
with open boundary conditions on the subsample and sam- 
ple. Note that the sample sizes are approximately separated 
by a factor of 2, except for the largest two sizes. For small 
h, the probability Pdo ~ 0.38, quantitatively consistent with 
a simple model with two states. For large h, Pdo — > as 
L —f oo, consistent with a single state. The solid line shows 
the step function that would obtain for Pdo in the oo-volume 
(and large w) limit, if it were the case that \m\ = 1 in the 
ferromagnetic phase, taking he — 2.270. The data is con- 
sistent with the calculated Pdo values approaching this step 
function at larger sizes L, for h ^ he- Note that at h — he, 
Pdo ~ 0.368 ± 0.006 < 0.379 . . . for L = 97 and L = 129, 
w = 3, consistent with a constant or slowly decaying value of 
Pdo at the critical point. The dashed curves are spline fits 
to organize the data visually. 



XI. SUMMARY 

In this paper, we have presented numerical results for 
the ground states of 3D random field Ising magnets fo- 
cusing on the transition between the ordered and disor- 
dered phases. Our results allow us to conclude that the 



transition is second order, though the magnitude of the 
magnetization vanishes very slowly as the critical random 
field strength, he, is approached from below. In addition 
to the magnetization, we have studied the stiffness of the 
system and some of the geometrical aspects, in particu- 
lar the fractal properties of domain walls at the critical 
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point. In general, the results agree very well with a sea! 
ing picture of the transition introduced some time agi 
and extensions of it to the properties studied here. 

Some earlier authors have suggested that the behavior 
of the RFIM near to the ordering transition will be more 
complicated than this scenario, arguing for some kind of 
"replica symmetry breaking" . Although, as is so often 
the case, the meaning of this term in this context has 
not been made clear, if we take it to imply the existence 
of many infinite volume ground states, this would have 
testable consequences. Although a full test of the de- 
pendence of the ground states on sequences of boundary 
conditions that this would imply is beyond the scope of 
today's computers and algorithms, we have made some 
preliminary tests on the dependence on boundary con- 
ditions. In particular, we have studied the probability 
that the configuration in a fixed volume at the center of 
a sample can be induced to differ from both the fixed 
s = -1-1 and fixed s = — 1 boundary conditions by various 
other boundary conditions. With the range of bound- 
ary conditions we have tested, this probability vanishes 
in the expected manner as L — > cx). Indeed, the power 
law dependence of this probability on L and the scaling 
with h — he are consistent with the domain wall fractal 
dimension and correlation length exponents determined 
by other methods. Our results are thus consistent with 
a single disordered to ordered transition at he, with a 
unique state in the disordered phases and a pair of states 
("up" and "down") in the ordered phase. Recent simula- 
tions at finite [temperature in smaller systems by Sinova 
and Canright,E3 who used the spectrum of the spin-spin 
correlation matrix and P{q) distributions, also suggest a 
single transition. 

At non-zero temperature, the thermodynamic proper- 
ties of the phase transition are believed to be similar to 
those at zero temperature: the transition is governed by 
a zero temperature fixed point. But at positive tempera- 
ture, one can also consider dynamic effects; indeed, as has 
been known for a long time, these dominate both Monte 
Carlo simulations and experiments. As first pointed out 
by Griffiths,LJ random systems can have singularities — 
albeit very weak ones — in thermodynamic properties 
well before the transition is reached and this will be the 
case for the RFIM. These rare region effects are unobserv- 
able as far as equilibrium properties iHrGlassical systems, 
but do have dynamic consequences EIE3'E3 In the Grif- 
fiths region above the transition, the average dynamic 
autocorrelations will decay more slowly than exponen- 
tially because of the effects of anomalously ordered lo- 
cal regions. Perhaps this kind of rare-region effect, and 
the more interesting but related effects that occur as the 
transition is approached, are all that is meant by "replica 
symmetry breaking". If this is the case, then it would 
be nice if the proponents of these ideas would say so. 
If not, then it is incumbent upon them to come up with 
some testable predictions. If these can be tested by static 
ground state properties, the RFIM is as good as system 
as any on which to perform such tests as the system sizes 



that can be studied are quite impressive: comparable to 
the largest that can be studied by Monte Carlo simula- 
tions in pure systems. 
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APPENDIX A: ALGORITHM 
IMPLEMENTATION 

We briefly describe here the algorithm and code used, 
including modifications to the conventional RFIM max- 
flow problem; outline the verification procedure for the 
code; and briefly outline the statistics and error bar pro- 
cedure. 



1. Base code and modifications 

There is a now well-known mapping of the RFIM 
ground state problem to a min-cut/max-flow problem.ta 
This correspondence and the push-relabel algorithm for 
the max-fiow problem, including terms used here (such 
as layers and excesses) , is well described in reviews and 
texts, such as Refs. 18, pq,_and pfl. The implementa- 
tion of the Goldberg- TarjadEB max-flow algorithm that 



we started with was the ili-fPrf code in C written by 
Cherkassky and Goldberg,E£l which in general performs 
quite well for a number of graph topologies. 

We modified the code to be more compatible with the 
CH — h language and developed objects (including sam- 
ples, configuration subsets as windows, and random num- 
ber generators) to conveniently implement a variety of 
boundary conditions and analyses. One very simple ben- 
efit of an integrated code is that the graph input, which 
is quite costly when read as a text file, is greatly sped 
up. More importantly, the short main routine was easily 
modified to compute answers to a wide number of ques- 
tions. 

The most significant change to the core push-relabel 
code was a modification that allowed for positive and 
negative excesses. This modification was developed in 
collaboration with D. McNamara.EEl The central idea is 
the elimination of the source and sink nodes, which con- 
ventionally have links to the nodes of the graph repre- 
senting the spins Si, in favor of introducing nodes with 
a negative excess. The first step in the conventional al- 
gorithm pushes as much flow as possible from the source 
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onto the lattice nodes. This step is replaced in our code 
with an initialization where a positive excess hi is placed 
on each node for which hi > 0. The connections to the 
sink are substituted for by placing a negative excess on 
the nodes with hi < 0. The push-relabel algorithm then 
proceeds with the usual steps, with the nodes with pos- 
itive excess having their excess pushed and their heights 
relabeled, as appropriate. The negative excess nodes act 
as sinks for the positive excess, until such a nodes total 
excess becomes positive. Besides removing the links to 
the source and sink, the global relabeling step must be 
modified. Instead of carrying out a breadth first search 
from the sink node, the breadth first search instead starts 
from the nodes with negative excess. (If no negative ex- 
cess nodes remain, the algorithm terminates with flow 
equal to the sum of the positive hi.) The initial totals of 
the positive and the negative excesses are compared with 
the final totals: the decrease in the total positive flow, 
for example, gives the maximum flow through the graph. 
The spin configuration and magnetization is determined 
by counting the number of nodes that are in the maximal 
layer. 

The removal of the source and sink nodes reduces the 
amount of memory used by an amount l/(d-|- 1) relative 
to the conventional memory requirements and results in 
a slight speedup. For the largest lattice sizes studied 
(256'^), memory requirements were reduced at the cost 
of speed. If pointers and integers each require 4 bytes, 
the Cherkassky and Goldberg implementation requires 16 
bytes for each arc and 32 bytes for each node (counting 
the layers as part of the per- node requirements.) The use 
of pointers was retained for system sizes up to 128'^. For 
a regular lattice, however, the nodes at the end of each 
arc, sister arcs, and the list of arcs at each node can be 
recomputed whenever needed. For a cubic lattice, then, 
the number of bytes per node is reduced from (6 x 16 -t- 
32) = 128 bytes to (6 x 4-1-32) = 56 bytes. For L^ = 128^ 
samples, the running time increased by a factor of « 2.5, 
primarily due to the recomputation of the tail nodes of 
the arcs and the sister arcs. 

One modification for the 256^^ samples was made that 
is not strictly sound, in that the algorithm could con- 
ceivably fail. In order to save memory, a limit to the 
maximum number of layers was implemented. In the 
Cherkassky and Goldberg code, the number of layers al- 
located is given by the number of nodes in the graph. 
In practice, however, far fewer are needed. A check over 
1000 samples for 8 < L < 128 was carried out, for several 
values of h. The number of layers needed, K , appears to 
be largest for h w he- At this value, the sample mean of 
the maximum layer needed is about K ^ 2L. The dis- 
tribution over samples of the number of layers needed 
roughly scales with L, though the dependence of the 
width of the distribution could be L^ , with x near 1. In 
any case, the distribution drops off very quickly with K. 
The maximum number of layers needed over 1000 sam- 
ples scales roughly linearly with L, K"^'^^ « 7L, where 
the maximum is for periodic boundary conditions and 



over a range of h, 2.0 < h < 2.5, with a peak in K'^'^^{h) 
near he ~ 2.27 (though the peak is slightly above he for 
smaller samples.) The mean number of layers fits rela- 
tively well to a scaling collapse, with a maximum value 
scaling consistent with L^'^^ (or even Ly/lnL), scaling 
about he ~ 2.27 with i^ w 1.35. We set 5 x 10"* as the 
maximum number of layers for all sample sizes, which is 
nearly 200 x L for the largest samples studied. This num- 
ber of layers was easily sufficient for all samples studied. 
The amount of memory needed for cubic lattices is then 
48 4- 0{L~^) bytes per node. 



2. Verification 

The modifications made to the base code, while the- 
oretically sound (except for the limit on the number of 
layers), could inadvertently introduce errors, due to er- 
rors in coding. We therefore verified the code against the 
Cherkasskpraud Goldberg codes h_prf and hi_pr (version 
3.3) codeslijO and a selection of other codes that were 
not based on a push-relabel algorithm. This was done by 
having the production code write out the list of the hi. 
A small program then converted these hi into arcs and 
nodes into a graph description DIMACS format, using 
the conventional representation with a source and a sink. 

These graph descriptions were then used as input to 
h_prf, hi_pr, and other codes, such— as those devel- 
oped in the First DIMACS Challenge.B The flows and 
the magnetization from these available algorithms were 
then compared, sample by sample, with the produc- 
tion code. The precise comparison was done for a few 
tens of samples with ferromagnetic coupling strength 
J = 10, 10^, . . . , 10^, relative disorder strength h/J = 
1,2.3,3,10 and system sizes L = 4,8, ...,128. The pro- 
duction code and the other algorithms agreed in all cases, 
except when the flow exceeded 2^^, which was generally 
the maximum possible flow in the available algorithms. 
The production code used here does not have that restric- 
tion, as flow computation is done at the end by comparing 
the initial and final positive and negative excesses, which 
were summed up as double precision quantities. The code 
was able to handle larger flows consistently, as could be 
verified by scaling h and J to large values (multiplying 
J and the hi by a factor of, say, 10^, and checking that 
the maximum flow increases proportionately and that the 
configuration is unchanged.) 

For efficiency, we have used an integer algorithm, with 
a resolution of 10^ by replacing the hi by random integers 
found by rounding Zi x 10** towards zero, with Zi Gaussian 
random variables with zero mean and unit variance and 
the exchange by J/h x 10^. We checked that reducing this 
resolution by a factor of 10 for selected measurements did 
not affect the computed averages, such as the stiffness in 
the isotropic and anisotropic samples for systems up to 
128'^. For a resolution of 10^, there were discrepancies 
outside of statistical errors, but these discrepancies could 
be consistently explained by the effects of rounding to an 
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FIG. 31: Elapsed time for computing ground states in the 
RFIM, plotted vs. h, for linear sizes L = 8, 16, . . . , 128. The 
"fast" algorithm is applied, with the larger memory require- 
ments, on a 766 MHz Pentium III processor. The peak times 
per spin scales nearly linearly with L. 



integer, which shifts the width of the distribution of hi 
by a small, predictable amount. 

We also tested the assumption that the cut found was 
unique (that is, that the ground state was unique, for a 
given sample.) Some accidental degeneracies were found, 
at the level of a fraction of spins « 0.9 x 10~^, for h near 
2, including he- This would result in the magnetizations 
being in error at the level of < 10~^, well within the 
statistical errors. Increasing the resolution by a factor of 
10 increases the running time by about 7% and reduces 
the fraction of degenerate spins to « 2 x 10"''. As the 
degeneracies were for the most part attributable to single 
spins, were rare, and did not affect any of the sample 
averaged results in the cases we tested, the integer scale 
of 10^ was more than sufficient for this study. 

We have verified that our choice of random number 
generator does not affect the results. Specifically, we used 
two generators for the computa tions of the magnetization 
and domains (defined in Sec. VIII) in 256^ samples at 
h — 2.27 w he and found the results to agree within 



statistical errors (the results reported pool the results 
from these generators.) We also checked the results from 
the two generators against each other for a larger number 
of smaller systems. 

Though quantities computed and the details of our 
interpretation differ from previous work, the numerical 
values for the same sample sizes and measurements (for 
example, magnetization and the-iaegest cluster sizes) are 
consistent with published data.llile^ 



3. Timing 

Consistent with similar optimization problems related 
to physical problems, the typical CPU time needed to 



find the ground state scales roughly as N^-^ near /ic. 
Roughly, it takes about 16-20 times longer to find the 
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FIG. 32: Statistics for the "operations", pushes (a) and 
relabels (b), performed in computing ground states in the 
RFIM, plotted vs. h, for linear sizes L = 8, 16, . . . , 128. The 
peak number of operations per spin scales nearly linearly with 
L (near he). At both high and low h, the number of these 
rearrangements scales nearly as L"'^. 



ground state each time the sample size is doubled, for 
hfn he. Using CC on a 400 MHz Sun UltraSparc II (the 
San Diego Supercomputing Center Sun HDSCIOOOO), a 
256'^ lattice required 913 MB of memory total for the 
graph data, the instructions, and the data structures re- 
quired for analysis. Running time for this size and this 
architecture averaged 1.8 hours per sample, for h — 2.27. 
Run times, normalized to the elapsed time per spin, for 
the larger memory algorithm, with the full data struc- 
ture, are plotted in Fig. ^. Clearly, the shape of the 
elapsed time vs. h sharpens some as L increases. The 
peak running time scales as ~ L^° over the scales L = 8 
to L = 128. Further details of the scaling of the run- 
ning time and connections between the algorithm and the 
physical concepts of ground state degeneracy and corre- 
lation length are described in Ref.pSl 
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